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Growing  applications  of  shape-memory  alloys  (SMAs),  which  possess  the  capabil- 
ity to  remember  their  original  shape  after  large  deformations  without  undergoing  perma- 
nent plastic  deformation,  demand  computationally  efficient  constitutive  models  and  the 
three-dimensional  finite-element  formulation  to  design  SMA-based  devices  in  the  me- 
chanics community. 

In  this  thesis,  the  finite -element  formulation  is  based  on  the  mixed  three-field  Fraeijs 
de  Veubeke-Hu-Washisu  (FHW)  variational  principle  to  treat  several  element-locking 
problems  such  as  shear  locking  in  bending  of  thin  structures;  and  volumetric  locking 
during  inelastic  flow  in  plastic  or  shape-memory  materials.  Because  of  the  volume- 
conserving  effects  of  shape-memory  alloys  under  transformation  flow  and  the  shear  lock- 
ing effects  in  thin  structures,  the  appropriate  finite-element  formulation  is  needed.  We 
used  the  solid-shell  finite  element  formulation  incorporating  the  enhanced  assumed  strain 
(EAS)  based  on  three-dimensional  three-field  mixed  formulation  to  take  account  of  the 
incompressibility  and  shear  locking  problems  in  the  element  with  appropriate  EAS  pa- 
rameters. The  element  formulation  has  only  a displacement  degree  of  freedom. 

The  constitutive  model  developed  is  thermodynamically  admissible,  i.e.,  the  model 


xi 


does  satisfy  Clausius-Duhem  inequality.  Based  on  the  thermodynamic  description,  a 
macroscopic  constitutive  relation  of  shape-memory  alloys  is  given  by  its  Helmholtz  free 
energy  which  is  decomposed  into  elastic  and  inelastic  parts,  in  analogy  to  the  conventional 
plasticity.  Transformation  function  is  also  determined  from  the  internal  entropy  produc- 
tion rate  and  transformation  flow  is  associated  with  transformation  function  in  stress- 
temperature  space.  The  direction  of  transformation  flow  rule  is  normal  to  the  transforma- 
tion surface  in  stress-space  so  that  the  transformation  surface  is  expanding  isotropically 
in  stress  deviator  during  the  martensite  transformation  deformation,  and  is  similar  to  the 
case  of  the  J2  flow  plasticity  of  isotropic  hardening  law.  In  the  discretized  constitutive 
model,  we  also  present  a return-mapping  algorithm  as  a special  case  of  the  closest-point 
projection  schemes  to  be  implemented  into  finite-element  formulation. 

Numerical  results  for  the  given  specimen  subject  to  the  tensile  are  also  provided 
to  show  qualitatively  the  ability  of  the  present  model  and  to  compare  quantitatively  with 
other  models  in  uniaxial  testing.  In  order  to  show  the  computational  efficiency  of  solid- 
shell  element,  we  test  our  constitutive  model  in  simply-supported  bending  problems  and 
in  thin  shell  structure  too.  In  all  the  tested  cases,  the  results  from  the  solid-shell  element 
are  compared  with  the  results  from  the  pure  displacement-based  twenty-node  element. 
For  the  given  tests,  the  proposed  model  shows  an  improvement  of  Newton-Raphson  it- 
eration schemes  to  solve  the  system  of  equations  and  results  in  numerical  efficiency  in 
computational  procedures. 
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CHAPTER  1 
OVERVIEW 


Over  the  last  few  decades,  shape-memory  alloys  (SMAs)  have  become  among  the 
most  attractive  smart  materials  in  research  and  applied  science  communities,  because  of 
their  unique  material  behaviors  and  broad  applications  in  such  diverse  areas  as  aerospace 
engineering,  automotive  industry,  and  biomedical  engineering.  Shape-memory  alloys  are 
a class  of  metal  compounds  that  possess  the  capability  to  remember  an  original  shape  even 
after  relatively  large  deformations  without  undergoing  permanent  plastic  deformation. 

The  discovery  of  the  shape-memory  effect  in  Ni-Ti  in  the  1960s  prompted  inten- 
sive research  activities,  leading  to  a variety  of  shape-memory  alloys  and  even  industrial 
products.  The  field  is  still  expanding  in  depth  as  well  as  in  width  and  new  theories,  new 
alloys,  and  new  applications  appear  frequently.  The  most  widely  used  shape-memory 
material  is  an  alloy  of  nickel-titanium  (Ni-Ti).  Nitinol,  one  of  the  most  well-known 
shape-memory  alloys,  is  capable  of  recovering  strains  in  the  order  of  10%  when  heated 
above  its  phase  transformation  temperature.1  Other  shape-memory  alloys  include  copper- 
aluminum-nickel  (Cu-Al-Ni),  copper-zinc-aluminum  (Cu-Zn-Al).  While  the  material  be- 
haviors in  shape-memory  alloys  have  been  investigated  based  on  most  of  the  experimental 
work  since  its  discovery,  when  it  was  attracted  by  engineering  and  science  researchers,  the 
computer  simulations  incorporating  the  nonlinear  constitutive  laws  of  shape-memory  al- 
loys have  only  been  traced  back  to  the  early  1990s.  Since  then,  many  constitutive  models 
have  been  developed  based  on  either  microscopic  or  macroscopic  viewpoint. 

1 Nitinol  standing  for  Nickel  Titanium  Naval  Ordnance  Laboratory,  was  discovered  to  possess  the 
unique  property  of  having  shape-memory  when  William  J.  Buehler,  a researcher  at  the  Naval  Ordnance 
Laboratory  in  White  Oak,  Maryland,  discovered  this  shape-memory  alloy  in  1961.  For  more  details,  see 
Buehler,  Gilfrich  and  Wiley  [1963], 
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1 . 1 Motivations  and  Objectives 

To  effectively  design  and  fabricate  shape-memory  alloys  devices  or  structures  re- 
lated to  shape-memory  materials  and  other  inelasticity,  one  must  evaluate  and  understand 
the  material’s  response  under  conditions  simulating  operating  environments.  In  other 
words,  modeling  of  constitutive  relations  is  of  significance  because  it  allows  designers  to 
quantitatively  and  qualitatively  predict  material  behaviors  and  then  produce  design  based 
on  the  calculated  numbers.  The  modeling  is  quite  difficult  since  the  main  features  involve 
nonlinear  mechanism  or  large  strains  of  material  behaviors.  Thus  the  usual  theory  of 
elasticity  is  not  applicable  and  one  must  turn  to  the  nonlinear  theory  of  complex  material 
responses  such  as  plasticity,  superelasticity,  viscoplasticity,  etc. 

From  the  continuum  mechanics  viewpoint,  a thermo-elastic  martensitic  phase  trans- 
formation in  the  material  is  responsible  for  the  extraordinary  properties  of  shape-memory 
materials.  Hence,  it  is  desirable  to  model,  based  on  the  theory  of  thermodynamically  con- 
sistent in  the  sense  of  its  compatibility  with  the  Clausius-Duhem  inequality  as  a special 
formulation  of  the  second  law  of  thermodynamics.  The  first  objective  is  to  develop  a 
3-D  general  constitutive  law  to  study  mechanical  behavior  in  the  range  of  superelasticity 
of  shape-memory  alloy  under  proportional  loading  conditions  for  shape-memory  alloys. 
Unfortunately,  a lot  of  attention  has  been  paid  to  metallurgical  research  of  shape-memory 
alloys  in  three  decades  and  less  attention  was  dedicated  to  shape-memory  constitutive 
modeling.  The  present  model,  however,  does  not  claim  to  be  a universal  model,  but 
only  one  contribution  to  modeling  of  shape-memory  effects  for  SMAs  and  to  lead  to  a 
guidance  of  modeling  a more  general  model  incorporating  plastic  effects  into  it.  The  fol- 
lowing properties  are  to  be  satisfied  in  modeling  of  constitutive  relations  and  to  improve 
the  computational  efficiency  of  finite-element  analysis. 

In  J2  type  flow  which  is  a function  of  only  the  second  invariant  of  stress  deviator, 
the  direction  of  transformation  flow  rule,  conserving  the  normality  rule  with  respect  to 


3 


the  transformation  surface,  is  to  be  the  stress  deviator  so  that  the  transformation  surface  is 
expanding  isotropically,  like  the  well-known  phenomena  in  classical  plasticity.  Influence 
of  the  second  and  third  stress  invariants  to  describe  asymmetry  in  tension  and  compression 
uniaxial  loading  also  needs  to  be  investigated.  The  discrete  constitutive  model  being 
implemented  into  finite-element  formulation  also  preserve  the  closest-point  projection 
schemes  to  obtain  the  consistent  tangent  modulus  in  algorithmic  computational  model 
so  that  it  will  show  faster  convergence  in  the  global  Newton-Raphson  nonlinear  iteration 
schemes,  compared  to  the  other  models. 

1 .2  Chapter  Overview 

Chapter  2 provides  a background  on  shape-memory  alloys  and  the  thermodynamic 
principles  because  the  present  model  is  to  be  thermodynamically  consistent  in  the  sense 
of  the  Clausius-Duhem  inequality. 

Given  the  variety  of  potential  usages  for  SMA,  the  ability  to  accurately  model  and 
analyze  via  finite  element  procedure  has  been  extremely  attractive.  The  incorporation  of 
SMA  finite  element  procedure  into  the  design  stages  of  new  products  could  reduce  devel- 
opment times  and  costs  dramatically.  Chapter  3 presents  the  finite  element  formulation 
of  the  solid-shell  element.  After  a presentation  of  the  kinematics  assumption  and  the 
formulation  based  on  the  Fraeijs  de  Veubeke-Hu-Washizu  (FHW)  variational  principle 
(Felippa  [2000]),  the  finite-element  discretization  together  with  the  enhanced  assumed 
strain  (EAS)  is  reviewed. 

Chapter  4 describes  a phenomenological  macroscopic  constitutive  model,  the  ther- 
modynamic consistent  in  the  sense  of  its  compatibility  with  the  Clausius-Duhem  inequal- 
ity, is  described  using  the  assumption  of  Helmholtz  free  energy  from  which  stress-strain 
relations  are  derived  and  the  dissipation  function  for  the  driving  force  of  phase  trans- 
formation. Significance  of  normality  condition  in  J2  associative  flow  rule  will  also  be 
described.  The  new  isotropic  transformation  hardening  function  is  shown  incorporating 
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linear  term  and  nonlinear  term,  the  Voce-type  saturation  strain  hardening,  which  is  widely 
used  in  classical  plasticity  (Voce  [1948]).  Influence  of  the  second  and  third  stress  invari- 
ants to  describe  asymmetry  in  tension  and  compression  uniaxial  loading  is  also  to  be 
investigated.  The  constitutive  model  being  implemented  into  finite-element  formulation 
is  also  to  be  preserved  the  closes-point  projection  schemes  to  obtain  the  efficient  con- 
sistent tangent  modulus  in  algorithmic  computational  model  so  that  it  will  show  faster 
convergence  in  the  global  Newton-Raphson  nonlinear  iteration  schemes. 

Numerical  examples  are  shown  for  the  several  loading  cases  undergoing  stress- 
induced  martensitic  phase  transformation  in  Chapter  5.  For  the  benchmark  test,  a simple 
bar  under  uniaxial  tensile  loading  case  is  tested  and  shows  the  numerical  efficiency  of 
the  present  constitutive  model.  In  three-point  simply-supported  bending  simulation,  the 
propagated  transformed  zone  is  shown,  and  it  is  similar  to  uniaxial  tensile  loading  case.  In 
order  to  show  the  computational  efficiency  of  solid-shell  element,  we  test  our  constitutive 
model  in  simply-supported  bending  problem  and  test  in  thin  shell  structure  as  well.  In  all 
the  tested  cases,  the  results  from  solid-shell  element  are  compared  with  the  results  from 
the  pure  displacement-based  twenty-node  element.  For  the  given  tests,  our  model  shows 
an  improvement  of  Newton-Raphson  iteration  schemes  to  solve  the  system  of  equations 
and  results  in  numerical  efficiency  in  computational  procedures.  Comparative  results  to 
other  models  are  shown  so  that  this  work  will  provide  valid  data  for  the  determination  of 
three-dimensional  constitutive  models  for  SMA  as  well  as  a framework  for  future  analysis 
of  SMA  through  numerical  procedures. 

Finally  Chapter  6 concludes  and  guides  the  work  with  future  directions  for  research. 


CHAPTER  2 

BACKGROUND  KNOWLEDGE 


2. 1 Typical  Features  of  Shape-Memory  Alloys 

Phenomenological  responses  in  shape-memory  alloys  (SMAs)  are  reasonably  well 
understood.  The  features  of  SMAs  can  be  described  by  variations  of  the  martensitic  frac- 
tion corresponding  to  the  content  of  martensite  during  martensitic  and  reverse  transfor- 
mation and  the  transformation  temperatures  shown  in  Figure  2.1.  If  one  takes  a stress-free 
SMA  at  high  temperature  in  the  austenite  (or  parent)  phase  and  cools  it,  a gradual  trans- 
formation to  the  low  temperature  or  martensite  phase  is  achieved.  Several  variants  of 
the  martensite,  including  multiple  twins  are  obtained  in  this  cooling  process.  Since  for 
a given  SMA,  the  transformation  from  austenite  to  martensite  (and  the  reverse  transfor- 
mation due  to  heating)  is  gradual,  an  important  state  variable  is  the  martensite  fraction 
(£)  with  0 < £ < 1 where  £ = 1 and  £ = 0 correspond  to  an  all  martensite  phase  and 
all  austenite  phase,  respectively.  Associated  with  this  state  variable  £ are  four  tempera- 
tures referred  to  as  martensite  finish  (My),  martensite  start  (Ms),  austenite  start  (As),  and 
austenite  finish  (Ay)  temperatures.  When  an  SMA  in  austenite  phase  (£  = 0)  is  cooled, 
Ms  is  the  temperature  at  which  the  transition  to  martensite  begins  and  My  is  that  tem- 
perature at  which  the  transition  is  completed,  i.e.,  £ = 1 so  that  the  material  is  composed 
completely  of  the  martensitic  variants.  An  analogous  description  defines  the  temperatures 
As  and  Ay.  Sometimes  it  is  possible  to  have  Ms  > As,  but  the  hysteresis  effect  is  small. 
2.1.1  Deformation  Mechanisms  of  Shape-Memory  Alloys 

Shape-memory  alloys  display  different  mechanical  behaviors  at  low  and  high  tem- 
peratures. These  unique  material  characteristics  can  be  described  by  crystallographically 
reversible  martensitic  transformation  between  the  phase  that  is  stable  at  low  temperature 
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Figure  2.1.  Variations  of  the  martensitic  fraction  during  forward  and  reverse  transforma- 
tions. 

called  martensite  and  the  phase  that  is  stable  at  high  temperature  called  austenite  in  the 
presence  of  external  applied  stress  and/or  changes  in  temperature. 

It  is  better  to  describe  the  way  the  atoms  are  arranged  in  general  materials.  The 
atoms  are  arranged  in  layers  and  rows  called  a lattice  (Figure  2.3).  Strong  forces  exist 
between  each  atom  and  its  neighbors.  These  forces  hold  the  atoms  together,  making  the 
material  very  strong.  When  the  material  is  stretched,  each  atom  is  pulled  slightly  away 
from  each  of  its  neighbors.  When  the  stretching  force  is  removed,  the  atoms  each  pull 
their  neighbors  back  into  place  (pseuo-plasticity  or  shape-memory  effects).  However,  if 
the  stretching  force  is  strong,  then  the  atoms  get  pulled  too  far  away  from  the  atoms  they 
were  originally  next  to.  When  the  stretching  force  is  taken  away,  the  atoms  stay  near  their 
new  neighbors,  instead  of  going  back  to  their  original  places  (plasticity ).  Of  course,  if  the 
stretching  force  is  big  enough,  the  atoms  will  be  ripped  apart,  and  the  material  will  break 
(failure).  See  Figure  2.2. 

In  Ni-Ti  systems  SMAs,  nickel  atoms  and  titanium  atoms  are  arranged  alternatively 
in  the  lattice  structure  in  the  parent  phase  state.  The  structure  of  arrangement  of  SMAs 
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Figure  2.2.  Stress-strain  curves  showing  PT  and  plastic  deformation. 

is  changed  because  the  layers  of  the  lattice  are  arranged  differently  dependent  on  temper- 
atures. Typically,  there  are  two  important  features  in  SMAs.  Figure  2.3  shows  that  the 
austenite  structure  (b)  can  be  converted  into  the  internally  twinned  martensite  structure 
by  cooling  or  supplying  a shear  stress.  If  the  temperature  is  below  the  austenite  transfor- 
mation initiation  temperature  (As),  SMAs  in  austenite  can  be  transformed  to  a martensite 
phase  whose  crystallo-graphic  structure  has  many  variants,  typically  sheared  platelets. 
Because  the  transformed  martensitic  phase  is  self-accommodating,  the  deformation  on 
transformation  to  martensite  is  zero.  The  martensite  deforms  by  a twinning  mechanism 
that  transforms  the  different  variants  to  the  variant  that  can  accommodate  the  maximum 
elongation  in  the  direction  of  the  applied  force.  The  interfaces  between  variants  in  the 
martensite  phase  slip  very  easily  and  the  material  is  deformed  at  low  applied  stresses. 
And  the  material  is  triggered  into  returning  to  its  original  shape  by  the  application  of  heat. 
That  is,  if  temperature  less  than  As,  the  reverse  transformation  will  not  take  place  after 
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unloading  and  the  deformation  will  be  kept  [(b)  — > (a)  —>  ( d ) — > (e)]  as  the  so-called 
residual  strain.  In  the  following  heating  process  (e)  — > (b)  in  Figure  2.3,  the  residual 
strain  will  disappear  and  the  original  state  is  restored  for  the  reverse  transformation.  The 
austenite  phase  has  only  one  possible  orientation,  thus  all  the  possible  deformed  struc- 
tures of  the  martensite  phase  must  revert  to  this  one  orientation  of  the  austenite  memory 
phase  and  the  material  recovers  its  original  shape.  It  is  called  the  shape-memory  effects 
(SME). 

SMA  also  display  superelasticity  (SE),  that  is  a mechanical  type  of  shape  mem- 
ory. This  effect  is  observed  when  alloys  are  considerably  deformed  at  temperature  above 
austenite  finish  ( A f).  When  the  stress  is  removed,  the  martensite  reverts  to  austenite  and 
the  material  recovers  its  original  shape.  This  process  (a)  — > (c)  in  Figure  2.3  shows 
the  martensite  transformation  induced  by  a cycle  of  stress  loading/unloading.  In  this  pro- 
cess, a significant  deformation  is  observed  for  the  uniformly  oriented  detwinned  deforma- 
tion and  this  process  is  an  important  feature  of  stress-induced  martensite  transformation 
in  SMAs.  That  is,  when  we  apply  a load  on  SMA  specimen  in  the  austenite  phase  at 
temperature  above  Af,  the  deformation  will  be  recovered  with  the  disappearance  of  the 
martensite  after  unloading  (c)  — > (a). 

2. 1 .2  Macroscopic  Phenomena  of  Shape-Memory  Alloys 

The  relation  between  stress  and  strain  in  this  loading/unloading  process  at  T > 
Af  is  shown  in  Figure  2.4-A  and  the  behavior  is  classified  as  superelastic  effect  (SE). 
Stresses  beyond  the  austenite  yield  limit  will  cause  detwinning  of  the  lattice,  forcing 
the  formation  of  martensite  directly  from  the  austenite  phase.  Thus,  the  SE  is  produced 
by  a stress-induced  martensite  forward  phase  transformation  (FPT)  and  its  reverse  phase 
transformation  (RPT)  process  (Figure  2.5). 

The  stress-strain  characteristics  of  austenite  Ni-Ti  and  martensite  Ni-Ti  are  shown 
in  Figure  2.4-B.  If  stress  is  applied  to  the  specimen  at  T < Mj,  when  the  alloy  is 
fully  martensite,  the  twinned  lattice  plate  will  first  deform  elastically  up  to  the  elastic 
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(a)  Martensite  A*)  Austenite  (c)  Martens^ 


(d)  Martensite  Martensite 


Figure  2.3.  Deformation  mechanism  and  lattice  structure  of  shape-memory  alloys 


limit  (martensite  yield  limit).  If  the  stress  exceeds  the  martensite  yield  limit,  detwin- 
ning  occurs  and  a large  non-elastic  deformation  entails  as  the  twin  boundaries  of  the 
lattice  shift  to  the  fully  detwinned  state.  Further  loading  causes  elastic  deformation  of 
the  detwinned  lattice.  Unloading  from  the  fully  martensite  state  will  result  in  residual 
deformation.  Unlike  in  plasticity,  the  shape-memory  materials  are  returning  to  the  orig- 
inal shape  by  application  of  heat,  shown  in  Figure  2.4-B.  It  is  called  the  shape-memory 
effects  (SME).  As  seen  in  the  above  simple  explanation  of  the  deformation  mechanism 
of  the  SMA,  the  deformation  behavior  or  the  shape-memory  effect  of  the  SMA  is  a com- 
bination of  the  temperature-induced  martensite  phase  transformation  (TIMPT)  and  the 
stress-induced  martensite  phase  transformations  (SIMPT)  between  the  austenite  structure 
and  the  martensite  structure. 

Since  SMA  have  unique  properties  which  ordinary  metals  do  not  have,  they  have 
high  potentiality  for  many  applications;  from  engineering  applications  and  biomedical 
implants  or  devices  to  everyday  consumer  products  (Duerig,  Melton,  Stokel  and  Wayman 
[1990],  Otsuka  and  Wayman  [1998],  Birman  [1997a]).  Pipe  couplings  are  the  first  most 
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Figure  2.4.  Stress-strain  curve  corresponding  to  stress-induced  martensite  phase  transfor- 
mation (A)  and  to  a combination  of  temperature-induced  martensite  transformation  and 
stress-induced  martensite  transformation  (B). 


Figure  2.5.  Stress-temperature  ( a — T)  phase  diagram  of  superelasticity. 


SMA  pipe  coupling 


✓ 


(a) 


Heating 


(b) 


Figure  2.6.  SMA  applications:  Pipe  coupling. 

successful  applications  of  SMA  to  hydraulic  systems  of  F-14  aircraft  in  1969.  In  pipe 
couplings,  SMA  tube  couplers  are  fabricated  being  cooled  to  cryogenic  temperatures  to 
produce  the  martensitic  phase  as  a ring  or  sleeve  with  the  inner  diameter  smaller  than  that 
of  the  pipes  to  be  connected.  After  expanding  at  an  ambient  temperature,  two  pipes  can 
easily  be  installed  into  the  pipe  coupling.  Once  the  coupler  joint  is  heated  to  T > As,  the 
pipe  coupling  begins  to  shrink  and  a surplus  fit  occurs  between  the  pipe  coupling  and  the 
pipes,  thus  producing  a very  effective  seal  together  (Fig.  2.6). 

Another  type  of  smart  materials  are  composites  of  SMA  embedded  in  non-SMA 
composite  structures  to  reinforce  and  to  be  used  for  vibration  control  of  the  composite 
structures,  since  the  elastic  constants  can  be  changed  widely  by  changing  the  temperature 
in  the  transformation  temperature  range.  The  SMA  fibers  or  wires  are  stress  loaded  and 
deformed  at  low  temperatures  (T  < As ) in  the  martensitic  phase.  They  are  then  unloaded 
to  generate  some  martensitic  residual  strain.  Once  embedded  in  the  elastic  matrix  they 
can  be  heated  using  an  electrical  or  thermal  input  to  achieve  residual  strain  recovery.  This 
recovery  can  be  free,  fully  restrained,  or  controlled  (Birman  [1997a]). 

Nitinol  is  used  in  a vast  variety  of  applications  in  medical  implants  and  instruments 
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because  of  its  excellent  corrosion  resistance,  mechanical  properties,  and  biocompatibility 
(Gil  and  Planell  [1999]).  Thus,  Nitinol  has  become  good  material  choice  for  desingers 
of  self-expanding  stents  and  orthodontics.  Self-expanding  stents  are  manufactured  with 
a diameter  larger  than  that  of  target  vessel,  crimped  and  restrained  in  a delivery  system, 
then  elastically  released  into  the  target  vessel  to  reinforce  blood  vessels  (Duerig,  Tolomeo 
and  Wholey  [2000],  Stoeckel,  Bonsignore  and  Duda  [2002]). 

In  mechanical  or  civil  damping  devices,  by  using  the  ability  of  SMA  to  absorb 
large  amounts  of  strain  energy  and  release  it  as  the  applied  strain  is  removed,  it  has  been 
used  for  damping  devices,  especially  in  seismic  structures  (Saadat,  Salichs,  Noori,  Hou, 
Davoodi,  Bar-on,  Suzuki  and  Masuda  [2002],  Dolce  and  Cardone  [2001],  DesRoches  and 
Delemont  [2002]). 

2.2  Principles  of  Thermodynamics 

Thermodynamics  deals  with  the  exchange  of  energy  in  its  various  forms,  thermal, 
chemical,  and  mechanical,  when  temperature  changes  are  significant.  In  the  sense  of  ma- 
terial behaviors,  thermodynamics  is  the  science  of  choice  for  constitutive  modeling  since 
it  limits  the  possibilities  for  material  behavior,  thereby  simplifying  the  specification  of 
constitutive  models.  These  restrictions  are  especially  useful  in  predicting  the  complicated 
thermo-mechanical  response  inherent  to  shape-memory  materials.  Based  on  the  thermo- 
dynamic description,  a macroscopic  constitutive  relation  of  shape-memory  alloys  is  deter- 
mined by  its  Helmholtz  free  energy  which  is  decomposed  into  elastic  and  inelastic  parts, 
in  analogy  to  the  classical  plasticity.  Transformation  function  is  also  determined  from  the 
internal  entropy  production  rate  and  transformation  flow  is  associated  with  transformation 
function  in  stress-temperature  space.  The  direction  of  transformation  flow  rule  is  normal 
to  the  transformation  surface  in  stress-space  so  that  the  transformation  surface  is  expand- 
ing isotropically  in  stress  deviator  during  the  martensite  transformation  deformation,  and 
it  is  similar  to  the  case  of  the  J2  flow  plasticity  of  isotropic  hardening  law. 
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The  deformation  in  shape-memory  materials  is  also  assumed  to  be  rate-independent 
because  the  temperatures  are  too  low  to  cause  significant  thermal  activation  (Boyd  and 
Lagoudas  [1996,  p.8 1 5]).  We  shall  also  suppose  that  the  process  of  deformation  of  a 
material  occurs  so  slowly  that  thermodynamic  equilibrium  is  established  in  the  body  at 
every  instant,  in  accordance  with  the  external  conditions. 

2.2.1  First  Law 

According  to  the  first  law  of  thermodynamics,  the  total  time  rate  of  change  in  inter- 
nal energy  and  kinetic  energy  balances  the  total  supply  of  energy  through  external  forces 
and  heat  (Maugin  [1992,  p.34],  Malvern  [1969,  p.226]).  In  other  words,  the  change  of 
kinetic  energy  and  internal  energy  in  a control  volume  is  equal  to  the  work  done  on  the 
surface  by  the  traction  and  by  the  body  forces  and  the  heat  convected  into  the  volume  and 
heat  created  by  body.  Let’s  consider  a domain  Q with  a boundary  d£l. 


where  the  internal  energy  per  unit  mass  (or  specific  internal  energy)  is  denoted  by  £ and 
its  kinetic  energy  is 


£ + K — Winput  + Q 


(2.1) 


The  internal  energy  £ is 


(2.2) 


(2.3) 


The  actual  power  input  of  external  forces  acting  on  the  body  is  written 


(2.4) 


where  b is  the  body  force  per  unit  mass,  t is  a traction  vector  per  unit  area,  and  v denotes 
the  velocity  field.  With  the  property  of  tt  — on  d(d£l)  and  the  identity  of  diver- 
gence theorem  by  which  the  surface  integral  may  be  transformed  to  the  volume  integral 
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of  d(ajiVi)/dxj,  the  power  input  of  external  forces  can  also  be  written  as: 

winput  = +Pbi)  + °j*^L]  (2.5) 

where  the  expression  in  parentheses  equals  p^fr  by  the  equations  of  motion.  Hence  for 
small  deformations  case 

W 'inpUt  = ~ j ^vv  dQ  + J a : e d.9.  (2.6) 

W^int 

where  Wint  is  defined  as  the  internal  work. 

Let  Q be  the  amount  of  rate  of  heat  supply  into  the  body;  it  consists  of  two  terms, 
the  heat  generated  within  the  volume  Q and  the  heat  received  by  conduction  through  the 
boundary  d£l  of  Q: 

Q=  [ Q dfl  = [ pr  dn  — [ q-nd(dfi)  (2.7) 

in  in  J act 

where  r is  the  internal  heat  source  per  unit  time  per  unit  mass,  q is  the  heat  flux  vector 
per  unit  area  through  the  surface  of  the  body,  and  n is  the  outward  unit  normal  to  dfl.  By 
using  the  divergence  theorem, 

[ q*ni(ffl)  = [ div(q)  dfl  (2.8) 

J dQ  J 

we  obtain  the  local  energy  balance  equation : 

p£  = a : e + Q (2.9) 

where  the  heat  received  by  conduction  through  the  boundary  of  Q is 

Q = pr  — div(q)  (2.10) 

In  component  form,  the  energy  balance  equation  can  also  be  written  as 

c • , dch 

p£  = + pr  - 


(2.11) 
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2.2.2  Second  Law 

In  addition  to  the  internal  energy  and  the  rate  of  heating,  it  is  necessary  to  introduce 
two  more  new  variables:  absolute  temperature(T)  and  entropy^).  Entropy  expresses  a 
variation  of  energy  associated  with  a variation  in  the  temperature.  The  second  law  states 
that  the  total  time  rate  of  change  in  entropy  is  always  greater  than  or  equal  to  the  rate  of 
supply  of  entropy  through  heat  (Maugin  [1992,  p.37]).  That  is,  mathematically 

V>  [ \ prdQ  - [ ^q-n  d{OQ)  (2.12) 

J n 1 J an  T 

The  above  relation  is  called  the  Clausius-Duhem  inequality.  The  total  entropy  is  defined 
as 

fj:=  / p$d«  = -Q  (2.13) 

if  p is  the  specific  entropy.  The  local  form  of  the  second  law  of  inequality  of  (2.12)  can 
be  written  as: 

PItj  PV  ~ ^Pr  + div(^)  > 0 (2.14) 

where  yv  is  called  the  internal  entropy  production  rate  per  unit  mass  so  that  fr/rj  can  be 
defined  as  the  internal  entropy  production  rate  per  unit  volume.  It  is  often  to  express  the 
entropy  by  dividing  the  reversible  part  and  irreversible  part: 

dp  = dp{r)  + dq[ir)  (2.15) 

where  ah/' ) denotes  the  part  of  this  increase  due  to  interaction  with  the  surroundings  of  the 
system  and  di y(jr)  denotes  the  part  of  this  increase  due  to  changes  taking  place  inside  the 
system,  i.e.,  the  inequality  implies  internal  entropy  production  in  an  irreversible  process. 
Then,  if  dp  denotes  the  heat  absorbed  by  the  system  from  its  surroundings,  we  have 

dv{r)  = Y-  (2.16) 

The  change  dp (tr)  is  never  negative.  If  dpd^  is  zero,  the  process  is  said  to  be  reversible. 
If  dpW  is  positive,  the  process  is  said  to  be  irreversible. 
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2.2.3  Clausius-Duhem  Inequality 

The  fundamental  inequality  containing  the  second  law  is  obtained  by  replacing  pr 
with  the  expression  resulting  from  the  equation  of  conservation  of  energy,  i.e.,  the  first 
law  of  thermodynamics  as  follows: 

PV  - j,[p£  - cr  : e + div(q)]  + div(^)  > 0 (2.17) 


Using  the  following  identity 


q d[qi/T]  1 dqi  dT  1 1 1 

d,V<-)  = = T + * &T  = Tdiv(q)  ~ T5q'Srad(T) 


(2.18) 


' T dxi 

and  multiplying  by  T we  obtain  the  new  form  of  Clausius-Duhem  Inequality  as  follows: 


p[Tq  - S\  + cr  : e - ^grad(T)*q  > 0 


(2.19) 


We  can  re-arrange  the  above  equation  as  follows: 

ip  '■=  Tpv  - p£  + cr  : e - ^grad(T)*q  > 0 (2.20) 

where  '{ p is  the  rate  of  internal  entropy  production  per  unit  volume  per  unit  temperature. 
The  term  ip  is  called  also  dissipation  rate  and  hence  is  equal  to  the  dissipation  power. 
It  is  determined  by  the  state  and  its  instantaneous  change,  i.e.,  by  the  state  variables  and 
their  derivatives.  We  thus  have 


ip  = V(afafT,T)>0 

where  V is  called  the  dissipation  function  and  ct " is  the  tensor  for  the  internal  variables 
(Ziegler  [1983,  p.62],  Malvern  [1969,  p.268]). 


CHAPTER  3 

FINITE  ELEMENT  FORMULATION 
3.1  Introduction 

In  the  last  several  decades,  the  development  of  advanced  analysis  techniques  for 
general  composite  and  sandwich  shells  has  been  of  interest  for  analyses  of  composite 
structures  covering  material  and  geometric  nonlinearities.  However,  the  computational 
implementations  of  suitable  shell  elements  have  continued  to  challenge  finite  element  re- 
searchers and  have  not  finished  yet  up  to  now,  especially  in  the  field  of  multilayer  compos- 
ite shells  and  non-linear  mechanics.  For  general  finite  element  analysis,  the  displacements 
are  interpolated  in  terms  of  nodal  point  displacements  in  the  formulation  of  the  contin- 
uum elements,  while  in  the  formulation  of  structural  elements  for  beam/plate/shell,  the 
displacements  are  interpolated  in  terms  of  midsurface  displacements  and  rotations.  Shell 
element  formulations  have  been  mainly  developed  within  the  context  of  the  so-called  de- 
generated shell  concept,  incompatible  modes  approach,  enhanced  strain  formulations,  and 
higher-order  theories  for  composites,  etc  (MacNeal  and  Harder  [1985],  MacNeal  [1998], 
Yang,  Saigal,  Masud  and  Kapania  [2000]).  Recently,  the  solid-shell  element  has  been  de- 
veloped by  Vu-Quoc  and  Tan  [2003]  and  has  been  well  proven  in  the  application  of  high 
element  aspect  ratio  without  locking  difficulties  of  element.  The  displacement-based  el- 
ement, however,  yields  very  poor  performance  in  the  bending-dominated  situations  and 
thin  shells.  The  solid-shell  element  has  the  same  displacement  dofs  as  in  regular  lin- 
ear eight-node  brick  solid  element,  which  avoids  the  need  for  introducing  finite  rotations 
as  variables  in  shell  problem,  and  does  not  require  the  use  of  transition  elements  when 
connecting  the  standard  solid  elements. 

In  order  to  obtain  the  same  performance  as  the  shell  element,  the  enhanced  assumed 
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strain  (EAS)  method  and  the  assumed  natural  strain  (ANS)  method  are  applied  in  the 
present  formulation.  To  improve  the  bending  behavior  of  low-order  elements,  the  EAS 
method  based  on  the  Hu-Washizu  functional  was  proposed  by  Simo  and  Rifai  [1990]. 
For  large  deformation  analyses,  the  Green-Lagrangian  strain  tensor  E and  its  conjugate 
Second-Piolar  Kirchhoff  stress  tensor  is  used  (Bischoff  and  Ramm  [1997],  Klinkel  and 
Wagner  [1997]).  The  ANS  method  is  also  employed  to  eliminate  the  locking  effects  from 
the  compatible  low-order  interpolations.  ANS  interpolation  is  the  most  successful  tool  to 
overcome  the  shear-locking  effect  in  the  4-node  displacement-based  shell  elements  (Mac- 
Neal  [1978],  Dvorkin  and  Bathe  [1984]).  We  apply  an  ANS  interpolation  of  the  compat- 
ible transverse  shear  strains  to  treat  shear  locking.  In  the  case  of  curved  structure  with 
geometric  nonlinearity,  there  is  another  locking  effect:  The  so-called  curvature-thickness 
locking  (Bischoff  and  Ramm  [1997]),  which  is  also  known  as  the  trapezoidal  locking  in 
piezoelectric  element  (Sze  and  Yao  [2000]);  this  type  of  locking  can  be  avoided  by  intro- 
ducing the  ANS  interpolation  of  the  compatible  transverse  normal  strain,  as  proposed  by 
Betsch  and  Stein  [1996]  and  Sprenger  and  Wagner  [1999]. 


In  this  section,  a brief  review  on  the  three-field  Fraeijs  de  Veubeke-Hu-Washizu 
(FHW)  variational  principle  (Felippa  [2000],  Fraeijs  de  Veubeke  [2001])  and  its  role  in 
the  EAS  method  will  be  discussed,  beginning  with  as  follows: 


where  the  displacement  u,  the  Green-Lagrange  strain  E,  and  the  second  Piola-Kirchhoff 
stress  S are  the  independent  variables.  In  (3.1),  the  expression  of  compatible  part  of 


3.2  Three-Field  Fraeijs-Hu-Washizu  Functional 


(3.1) 
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strain  Ec  (it)  in  terms  of  u is  given  in  (3.40)  and  (3.43);  Ws  is  the  stored  strain  energy, 
t the  traction  vector,  it*  and  t*  the  prescribed  displacement  on  the  boundary  Su  and  the 
prescribed  traction  on  the  boundary  Sc,  respectively.  All  variables  here  are  expressed 
in  the  initial  configuration  B0,  with  Sa  U Su  = dB0,  i.e,  the  closure  of  the  union  of  the 
boundary  Su  and  the  boundary  Sa  forms  the  complete  boundary  of  t30. 

The  next  step  in  the  EAS  method  is  to  introduce  an  enhancement  (or  enrichment) 
to  the  Green-Lagrange  Strain  Ec  and  of  the  variation  5EC  as  follows 


E = Ec  + E , 

(3.2) 

SE  = SEC  + SE  , 

(3.3) 

where  E and  SE  represent  the  enhanced  strain  and  its  variation,  whereas  E and  SE 
represent  the  enhancing  strain  and  its  variation.  Introduction  of  (3.2)  into  (3.1)  yields 

n (it,  E,s)  = f W8  ( Ec  {u)  + E)dV  - f S : EdV 

B0  Bo 

- ( u*  (b*  - u)  pdV  + f (u*  — u)  mt  ( S ) dS 

Bo  Su 

- j u‘t*dS  . (3.4) 

where  it  is  noted  that  the  traction  t is  actually  a function  of  S (Malvern  [1969,  p.69]),  and 
not  an  independent  variable. 

By  designing  the  approximation  for  the  stress  field  S and  the  approximation  for  the 
enhancing  strain  Ec  such  that  the  following  orthogonality  condition  holds: 

j S *.  EdV  = 0 . (3.5) 

B0 

then  the  number  of  independent  variables  in  the  functional  fl  in  (3.4)  is  reduced  to  two. 


n (u,  e)  = Jws  ( Ec  {u)  + e)  dV 

Bo 


that  is 
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Bo  Su 

leading  to  the  following  total  variation 


I u'b* pdV  + J (u*  - u)  •t(S)dS  - j u't*dS  , (3.6) 

£o  Su  S ^ 


511  (u,  E^j  = SUstiff  (u,  E^j  + 6Uext  (u)  = 0 , 


(3.7) 


where  from  (3.23)  and  (3.28),  we  have 

6Ustiff  (u,  E)  = f (5EC  (u)  + SE)  : ~ Ws  [Ec  (u)  +E)dV, 

B0 

SUext  = — J Sumb*pdV  — J Su*t*dS  . 

Bo  S„ 

Consider  a perturbation  of  the  displacement  u as  follows: 

ue  ( X , t)  = u (X,  t ) + eSu  ( X ) , 


(3.8) 

(3.9) 


(3.10) 


since  the  variation  Su  is  only  a function  of  space.  Therefore  the  variation  of  Id  in  (3.6) 
with  respect  to  the  displacement  u is 


d 

de 


n(U„B)j  = /^S  ^M.S„yV  - j Sn-b'pdV 

£=0  So  V / Bo 

— j SwtdS  — j 5u‘t*dS  , VSu  , (3.11) 


c.  J J dW8  . 

Since  the  second-order  tensor  — — is  symmetric,  the  first  term  in  (3.11)  can  now  be 

oh/ 

rewritten  as  follows 

dws,(dEc(u)  \ aws, 

'Su  ) = F . grad(dii) 


dE  \ du 


dE 


= Div  [FiiwSu) _ Div  (Fiw)  'Su  • <3-I2) 


where  the  first  equality  is  shown  in  Remark  3.1,  and  the  second  equality  is  obtained  with 


Leibniz  rule. 
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Remark  3.1.  With  the  deformation  gradient  F (u)  expressed  with  respect  to  the 
basis  ej  0 Ej  (Figure  3.2),  we  obtain  the  following  relations 


i rpj 


F = F)et  ®EJ  , F1  = F)E 


(3.13) 


F1  F = FlIgijFjEI  0 EJ  , with  gl3  = ei‘ej, 

) EJ6u  . 


SF  = . p, 

du  du 


(3.14) 

(3.15) 


Note  that,  if  {e^  e2,  e3}  were  a system  of  orthonormal  vectors,  we  would  have  gij  = Sij, 
however,  in  general,  it  is  not  necessary  to  have  {elt  e2,  e3}  orthonormal.  In  this  remark, 
we  retain  the  notation  g13  for  generality.  We  define  a symmetric  second-order  tensor  S as 


S :=  = SabEa  0 Eu  . 


dE 


(3.16) 


By  using  (3.16),  (3.40),  and  (3.14),  we  can  then  rewrite  the  first  term  in  (3.1 1)  as  follows 


dW • • dE°-5u  = s : 


<9£  du 


2 du 


' dF\ 


Ml 

y 


(5-u  . 


Since  from  the  symmetry  of  5,  (3.17)  becomes 

. r«_  rn  dF3B 

dE  * ^ ’ 


(3.17) 


(3.18) 


With  the  following  component  form  of  FS  using  (3. 1 3)i  and  (3.16),  that  is,  we  obtain 


FS  = F)S  ei  0 Eb  , 

and  then  the  following  expression 

r)F  f)F3 

FS  : SF  = FS  : Su  = F\SAB  gij^-Su  . 

au  J du 

Comparing  (3.20)  to  (3.18)  and  using  SF  = grad  (Su),  we  arrive  at 

dws 


(3.19) 


(3.20) 


dE 


: sec  = fs  : sf  = fs  : grad  (5 


u 


(3.21) 
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Using  the  divergence  theorem  and  Su  = 0 at  Su  , the  integration  of  the  first  term  in 
(3.12)  becomes 


/Div( 


.dWs 

dE 


= J ( F--J^Su  | • ndS 


(3.22) 


Bo  ' / dB0 

where  n is  the  outward  normal  vector  to  the  surface  dI30.  With  (3.12)  and  (3.22),  we  can 
write  the  variation  of  II  in  (3.1 1)  with  respect  to  u as 


de 


:n  (u„e)  = [ 

6 = 0 


n F 


.dWa 

dE 


dB0=SaL)Su 


-J 

Bo 


• SudS  — / Div  F 


.dWs 

dE 


•SudV 


- J 5u’b* pdV  - J 8u*tdS  - J Su‘t*dS 

Bo  Su  Sr 

"J 

Su 

■/ 


n*  F- 


Div  F 


.dWs 

DE 

dW„ 


•SudS  + / 

( ^dWs\ 
n*  F-—^-  - r 

J 

[ \ dEj  J 

■ SudS 


•5udV  = 0 , M5u  , 


(3.23) 


thus  the  Euler-Lagrange  equations  resulting  from  the  variation  of  II  with  respect  to  u in 
above  (3.23)  are 

.dWx 


t = n'\F- 


t*  = n-  F 


dE  y 
dW„ 


dE 


Div  ( F dFl 
\ dE 


on  Su  , 
on  , 
pb*  = in  Bo  , 


(3.24) 

(3.25) 

(3.26) 


Then  make  the  variation  of  II  of  (3.6)  with  respect  to  E by  the  perturbation 

Ee  = E + eSE  , 

it  follows  that 


(3.27) 


de 


n (uEs)  = J ( ^ - s ) : 5EdV  , VSE  , 

e=0  B0  ' ' 


(3.28) 
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the  Euler-Lagrange  equation  associated  with  (3.28)  is 

s _ 


dE  ’ 


(3.29) 


3.3  Kinematics  of  Solid-Shell  in  Curvilinear  Coordinates 


The  shell  kinematics  of  deformation  is  described  by  using  the  position  vectors  of 
a pair  of  material  points  at  the  top  and  at  the  bottom  of  the  shell  surface.  In  this  kine- 
matic description,  a straight  transverse  fiber  before  deformation  remains  straight  after 
deformation.  Such  transverse  fiber  need  not  be  normal  to  the  shell  mid-surface  before 
deformation,  as  well  as  after  deformation.  The  initial  (undeformed)  three-dimensional 
continuum  of  the  shell  geometry  (Figure  3.1)  is  described  by 


\ K1 + e3)  xu  + (i  - e)  xt  (e,  e)' 


i 


2 + Xl]  + [-^U  - Xl)  , 


£ = (e1,  £2,  £3)  e □ :=  [-1,  l]  X [-1,  l]  X [-1,  l] , 


(3.30) 


where  X (•)  is  the  mapping  from  the  biunit  cube  □,  parameterized  by  the  material  co- 
ordination (£\£2,£3),  to  the  initial  configuration.  The  image  X (£)  of  a point  £ = 
(£\£2,£3)  € □ is  next  represented  by  a linear  combination  of  the  position  vectors 
Xu  (*,  •)  and  X i ( * , * ) of  a point  in  the  upper  material  surface  £3  = +1,  and  a pint 
in  the  lower  material  surface  £3  = —1,  respectively. 

Similarly,  in  the  current(deformed)  configuration,  the  geometry  of  the  solid  shell  is 
described  by 


*(«  = \ h (?  ■.  e) + x,  (?  ■,  f)]  + \e  [*„  (e,  e)  - *,  (e,e) 


1 r 


2 K1  + £3)  xu  (e1,  £2)  + (1  - £3)  xt  (e,  £2)]  , 


£ = £1,£2,£3  €□, 


(3.31) 
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Figure  3.1.  Initial  (undeformed)  configuration  of  solid  shell:  convective  coordinates  f* 
and  position  vectors  Xu  and  X,. 

where  x (•)  is  the  mapping  from  the  biunit  cube  parameterized  by  (f1,  £2,  £3)  (material 
configuration),  to  the  current  (deformed)  configuration,  xu(-,  •)  and  xt  (•,  •)  the  posi- 
tion vectors  of  the  deformed  upper  and  lower  surface  of  the  solid  shell,  respectively. 

The  initial  configuration  is  related  to  the  deformed  configuration  (Figure  3.2)  by  the 
displacement  field  u as  follows 


The  convected  basis  vectors  G*,  i = 1,2,  3,  in  the  initial  configuration  B0  are  related  to 
the  position  vector  X and  the  converted  coordinates  £*  by 


*(0  = x(0  + ti(0. 


(3.32) 


(3.33) 


and  satisfy  the  following  relations 


(3.34) 
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X 


Figure  3.2.  Solid-shell:  Material  configuration  B,  Initial  configuration  B0,  and  Deformed 
configuration  Bt. 


where  Gy  are  the  components  of  the  metric  tensor  with  respect  to  the  basis  Gl  ® Gj 
in  the  initial  configuration  B0.  configuration.  To  simplify  the  presentation,  we  will  omit 
the  argument  (£)  and  simply  write  G*  and  G\  The  covectors  G 1 can  be  obtained  by  the 
following  relation 


Gl  = GijGj , with  [Glj]  = [Gy]"1  G R3*3  . (3.35) 


similarly,  the  convected  basis  vectors  gx  in  the  current  configuration  Bt  are  obtained  using 
(3.32)  and  (3.33)  as  follows 


dx  du 

~Wi=Gi  + W 


* = 1,2,3, 


(3.36) 
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and  satisfy  the  following  relations 


9i  '9J  = •fl'j  = 9ij  , *,  j = 1,  2, 3 , 


(3.37) 


where  are  the  components  of  the  metric  tensor  with  respect  to  the  basis  gl  <g>  g3  in  the 
current  configuration  Bt.  The  covectors  gl  can  be  obtained  as  follows 


9l  = 9lJ9j  , with  \gtJ  = [gl3]  € R 


.3x3 


(3.38) 


The  deformation  gradient  expressed  in  converted  basis  vectors  g{  and  Gl  takes  the 


form 


(3.39) 


Using  (3.39),  (3.34),  and  (3.37),  we  can  then  write  the  compatible  part  Green-Lagrange 
strain  tensor  with  respect  to  the  convective  coordinates  as  follows 

E°=\  (- eTe  = \ [{°l  ® 9i)  (?j  ® &)  - Gl3G'  ® Gj 

= \ (. 9ij  ~ Gn ) Gl  ® G3  = EtjG?  ® , (3.40) 

where  is  the  covariant  components  of  the  strain  tensor  Ec,  and  J2  the  identity  tensor 
expressed  with  respect  to  the  convected  basis  Gl  (8)  G 3 using  the  relations  ep  = GpGl 
and  eq  = GQjG 3 as  follows 

h = Spqep  ®eq  = 5pqGpG]Gl  ® Gj  = G^G1  ® G3  = GijGi  ® G 3 , (3.41) 


Using  (3.36),  the  metric-tensor  components  g^  in  (3.37)2  can  be  expressed  in  terms 
of  the  convected  basis  vector  Gt  and  the  displacement  vector  u as  follows 

ft,  = + + ^ 

du  du  du  du 

- G^G'-W  + WGi  + WW 


(3.42) 
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Substituting  (3.42)  into  (3.40),  the  covariant  components  of  the  compatible  Green- 
Lagrange  strain  tensor  Ec  read  as 


du  du 

9V+W'  l + 


du 


(3.43) 


The  second  Piola-Kirchhoff  stress  tensor  S is  conjugated  with  the  Green-Lagrange 
Strain  Ec,  can  be  expressed  with  respect  to  the  basis  G , <g>  G}  as  follows 


S = SijGi  ® Gj  , 

with  Sli  being  its  contravariant  components. 


(3.44) 


3.4  Finite  Element  Discretization 


In  this  section,  we  present  the  weak  form  and  finite-element  approximation  of  the 
proposed  solid-shell  element.  Let  the  initial  configuration  Bq  be  discretized  into  nonover- 
lapping  nel  elements  Bq}  with  numnp  nodes,  such  that  B0  « (J  Let  h denotes  the 

e=l 

characteristic  size  of  the  finite  element  discretization.  In  the  element  domain  B(0e\  the 
displacement  u,  its  variation  5u,  and  increment  Au  are  interpolated  as  follows 

U ss  uh  - N (£)  d(e)  , (3.45) 

8u  « 8uh  = N (£)  5d(e)  , Au  « Auh  = N ($)  Ad(e)  , (3.46) 

where  ATis  a matrix  containing  the  basis  functions  restricted  to  element  Bq\  and  d(e)  E 
r3x numnp  a matrjx  containing  the  nodal  displacements.  In  what  follows,  for  simplicity, 
we  will  omit  the  superscript  h in  uh,  and  simply  write  u. 

Within  a typical  element  (e),  the  variation  and  the  increment  of  the  compatible  part 
Green-Lagrange  Strain  Ec  is  related  to  the  variation  and  the  increment  of  displacement, 
respectively,  based  on  (3.43),  as  follows 


KL,  = B (dW)  • {^L,  = B (d<’>)  A d<*) 


(3.47) 
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where  the  components  of  Ef-  had  been  arranged  into  a 6 x 1 column  matrix  according  to 
the  Voigt  ordering 

= {-^Tl>  -^22  > 2-£'i2j  2£l23>  2£^i3}  , (3.48) 

and  where  B is  the  deformation-dependent  displacement-to-strain  operator,  where  de- 
tailed expression  is  given  in  the  Appendix  in  Vu-Quoc  and  Tan  [2003]. 

We  denote  the  admissible  variation  of  the  element  EAS  parameter  column  matrix 
a(e)  e Rneas  assocjateci  wjth  the  enhancing  strain  j by  5a(e),  where  neas  is  the 
number  of  EAS  parameters  Interelement  continuity  is  not  required  for  the  enhancing 
strain  E,  where  components  can  be  approximated  via  an  enhancing  strain  interpolation 
matrix  Q and  the  element  EAS  parameters  a. an  interpolation  applies  to  the  variation 
and  the  increment  of  E,  that  is 


= 0(£)«(e) 

(3.49) 

{A  Ea) 

l 1 6x1 

= £(OAc*(e>. 

(3.50) 

The  number  of  internal  parameters  a(e)  and  the  interpolation  matrix  Q will  be  discussed 
in  Section  3.6. 

The  consistent  tangent  operator  in  the  Newton  solution  procedure  is  constructed  by 
taking  the  directional  derivative  of  the  weak  form  at  a configuration  in  the  direction 
of  the  increment  A ^)U , where  the  left  subscript  k designates  the  iteration  number.  The 
tangent  operator  can  be  viewed  as  the  summation  of  the  material  and  geometric  tangent 
operators.  The  geometric  part  results  from  taking  the  variation  of  the  geometry  while 
holding  the  material  constant,  whereas  the  material  part  results  from  taking  the  variation 
of  the  material  while  holding  the  geometry  constant. 

Applying  a standard  finite-element  procedure  to  discretize  the  weak  form  (3.7),  we 
obtain  the  following  expression  at  the  element  level 

■5n<‘>  = sn%  + mS  = o , 


(3.51) 
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with  the  stiffness  part  SU^  for  the  static  case  (u  = 0 ) from  (3.8)  and  the  external-force 
part  SU^t  from  (3.9)  written  as  follows 


= j {S'S}dV+  j S{^}T  {S'’}dV  , (3.52) 

6oe>  B<e> 

fin®  — — J Swb*pdV  — J 8u‘t*dS  , (3.53) 

w(e)  o(e) 

°0  ‘-’a 

where  we  used  the  symbol  definition,  which  corresponds  to  the  second  Piola-Kirchhoff 
stress. 


S := 


dWs 
OE  ’ 


(3.54) 


to  alleviate  the  notation,  and  where  the  column  matrix  {5^}  has  its  coefficients  arranged 
in  the  same  Voigt  ordering  as  in  (3.48) 


{sii}  = [s11,s22,s12,s33,s23,s13' 


(3.55) 


It  should  be  noted  that  the  symbol  S in  (3.54)  is  simply  used  in  the  place  of 


dWs 

dE 


Remark  3.2.  The  linearization  of  the  weak  form  <511  (it,  E^j  can  be  accomplished 
by  the  truncated  Taylor  series  about  the  k — th  iterate  ( , ^E  j : 

8U  ((fc+1)u,  (fc+i )E)  ~ SU  ((fc)ii,  (k)E^ 


+ 


d(sn) 


d (it,  E^j 


(w=  (k)U , E= 


(Au.  A e) 


= ( {k)u  , {k)E  ) + V (AIT)  ( {k)u , {k)E  ) • (a u,  A e)  , (3.56) 


where  A u = (k+\)U  — ( k)U , A E = (k+\)E  — (k)E . To  compute  the  increment 
^Ait,  A E^j  in  the  Newton’s  solution  process,  we  simply  set  the  expression  in  (3.56)  to 


zero. 
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In  order  to  alleviate  the  notation,  we  will  omit  the  left  subscript  k designating  the  iterative 
index.  Using  the  approximation  (3.45),  (3.46),  (3.47)  and  (3.50)  in  (3.56),  the  increments 
A and  Aa^  can  be  computed  in  the  Newton’s  solution  process,  as  mentioned  above, 
by  using  the  following  equation 

v(6U{eA  (d(e),a^)  • (Ad(e\Aa^)  = + . (Ad(e)  Aa(e)\ 

V M J V ) 3(dMaw)  V J 


= . (Ad(e)  Aot(e)\  = _ f8 

d(d 1 j ^ 


H%  + 6U^t)  .V^.WaM,  (3.57) 


in  which  the  variation  in 

(3.52)  and  SU^t  in  (3.53)  now  take  the  form 

(d<‘> , «<'>)  = 

S<t‘)Tf%r  + Sa MTf™s , 

(3.58) 

II 

— fe 

*K 

-a 

f Br  [s13)dV  , 

fEAS=  f 9T{Sii}dV, 

(3.59) 

& 

<> 

inS(dW)  = 

-6dWf%  , 

(3.60) 

with  = 

J NTb*  pdV  + J 

NTt*dS 

(3.61) 

w(e)  o(e) 

°0  &<r 

Thus  the  left  hand  side  of  (3.57)  becomes 

9 (m%)  , 


da  (e) 


Ac*(e) 


17  («m^)  • (A d^\  AaW)  = ^^^A^ 

= [Sd{e)Tk^  + Sa^Tk^\  -Ad^  + [dd(e)Tk^a  + SatWkM] 


= 6d{e)T  [feg Ad(c)  + fcie>Aa<e)]  + Sa^T  [k& Ad(e)  + k^Aa® 

Let  the  matrix  of  tangent  elastic  moduli  D be  defined  as 

'dSVJ 


D 


jjijkl 


dE> 


kl 


G R 


6x6 


(3.62) 


(3.63) 


where  DljkL  are  the  components  of  constitutive  tensor  V in  the  convected  basis,  and  are 
subsequently  arranged  in  the  matrix  D according  to  the  ordering  of  the  strain  components 
in  (3.48)  and  of  the  stress  components  in  (3.55). 
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Using  (3.62),  (3.59)i,  (3.63),  and  (3.47),  we  obtain  the  following  expressions 


= Agf  = / (gt$  + btdb)  dv  , 


<9d(e) 


*2  = Ilf  = / , 

R<«> 


3(e) 


(3.64) 


(3.65) 


where  the  matrix 


G := 


SB  (d(e)) 
dd{e) 


(3.66) 


which  is  a function  of  the  coordinates  £,  and  the  stress  matrix  S which  is  related  to  the 
matrix  {SlJ}  in  (3.55),  given  in  the  Appendix  in  Vu-Quoc  and  Tan  [2003],  From  (3.62), 
(3.59)2,  (3.63),  (3.66),  and  (3.47),  we  obtain  the  remaining  parts  of  the  stiffness  matrix 


*£!  = =[fcS]T=  / 5rDBdV, 


»(e) 


I,(e)  _ 


5/ 


(e) 

£45 


= / QrDQdV  . 

„(=) 


(3.67) 


(3.68) 


It  follows  from  (3.57),  (3.58),  (3.60),  and  (3.62)  that  the  discrete  linearized  system 
of  equations  to  solve  for  the  increments  A and  is  given  by 


fcSAd<')  + fcWAaW  = 

k£>Ad(e)  + fc«Aaw  = 


/(e) 
7 ext 


/(e) 

7 stiff  1 


/(e) 

"7  £45  » 


(3.69) 

(3.70) 


or  in  matrix  form  as 


V.{e)  iu(e) 

™UU  ™UOl 

kie}  klei 


A d[e) 

Aa^ 


/(e)  _ /(e) 

7 ei«  7 siiff 
/(e) 

7 £45 


(3.71) 


Since  the  enhancing  strain  E is  chosen  to  be  discontinuous  across  the  element 
boundaries,  it  is  possible  to  eliminate  the  EAS  parameter  increment  Aa(e)  at  the  element 
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level,  before  proceeding  to  assemble  the  element  matrices  into  global  matrices.  Solving 
for  the  increment  Aa*e*  using  (3.70) 

A«w  = - [*£]'*  (f(As  + *£Ad<«>)  , (3.72) 


then  substituting  (3.72)  into  (3.69),  we  obtain  the  following  condensed  symmetric  element 
stiffness  matrix  k^  and  the  element  residual  force  vector  r ^ 


k[e) 

r{e) 


An  assembly  of  the  element  matrices  and  r ^ leads  to  the  global  system 


(3.73) 

(3.74) 


KTAd  = R,  (3.75) 

Tld  7ld 

with  Kt  = A.k$  , R = A , (3.76) 

e=l  e=l 

where  A denotes  the  finite-element  assembly  operator.  The  incremental  displacement 
Ad  can  be  solved  by  using  (3.75),  and  the  displacement  d and  d(e)  updated.  With  (3.72), 
the  incremental  displacement  A d^  is  used  to  compute  the  increment  Aa(e\  which  is  in 
turn  used  to  update  the  EAS  parameter 

3.5  Transverse  Shear  and  Normal  Strain  Treatment 

To  avoid  shear  locking,  we  use  the  assumed  natural  strain  (ANS)  method  by  Dvorkin 
and  Bathe  [1984]  for  the  transverse  shear  strains.  To  remedy  the  curvature-thickness  lock- 
ing (Bischoff  and  Ramm  [1997]),  we  adopt  the  ANS  method  by  Betsch  and  Stein  [1996] 
for  the  transverse  normal  strain.  The  assumed  natural  strain  (ANS)  method  was  originally 
prepared  to  relieve  the  shear  locking  problem  that  typically  arises  as  the  thickness  of  the 
shell  goes  to  zero  ( Hughes  and  Tezduyar  [1981],  Dvorkin  and  Bathe  [1984]),  and  was 
later  given  a mixed  variational  foundation  (Simo  and  Hughes  [1986]). 
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To  avoid  shear  locking,  we  adopted  the  ANS  method  as  applied  to  the  four-noded 
shell  element  in  Dvorkin  and  Bathe  [1984],  Here,  a linear  interpolation  of  the  compatible 
transverse  shear  strains  E[3  and  E £3  in  (3.40),  evaluated  at  the  four  midpoints  A , B,  C , D 
of  the  element  edges,  at  £3  = 0 (Figure  3.3),  is  applied  (Betsch  and  Stein  [1996],  Klinkel, 
Gruttmann  and  Wagner  [1999]) 

EtiNS  If  (1  - ?2)  Ef3  (£.„)  + (1  + a Bf3  (tc)  1 

E&NS  I \ (1 -?')£&  («D)  + (l+f‘)  / ’ l'  ’ 

where  the  coordinates  of  points  A,  B , C,  D are  £A  = (0,  — 1,  0),  = (1,  0,  0),  £c  = 

(0, 1,  0),  = (—1, 0, 0),  respectively. 

In  the  case  of  curved  thin  shell  structures  or  in  the  geometric  nonlinear  analysis, 
in  order  to  circumvent  the  locking  effect  from  parasitic  transverse  normal  strain,  we  em- 
ploy an  assumed-strain  approximation  for  the  covariant  component  E33  of  the  compatible 
Green-Lagrangian  strain  tensor,  as  done  in  Betsch  and  Stein  [1996]  for  a stress-resultant 
shell  formulation  and  for  eight-node  brick  element  in  Sprenger,  Gruttmann  and  Wagner 
[2000].  Here,  a bilinear  interpolation  of  the  transverse  normal  strains  sampled  at  the  four 
comers  E , F,  G , H of  the  element  midsurface  (Figure  3.3)  is  imposed,  that  is 

E^s  = E*i(e‘  ?2)  £3e3  («.)  . (3.78) 

i- 1 

with  TVj  =2(1+  £*  £J)  (1  + £2),  and  the  coordinates  of  the  corner  points  E , F,  G,  H 

being  = iE  = (-1, -1,0),  £2  = £F  = (1, -1,0),  £3  = = (1,1,0),  ^4  = 

*H  =(-1,1,0). 

3.6  Interpolation  of  the  Enhancing  Strain 

In  this  section,  we  first  review  the  regular  enhanced-strain  method  (Klinkel  et  al. 
[1999])  and  establish  the  optimal  number  of  internal  parameters  for  the  enhancing  strains 
in  the  present  solid-shell  element  to  pass  the  membrane  patch  test  and  the  out-of-plane 
bending  patch  test  (MacNeal  and  Harder  [1985]). 
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Figure  3.3.  Eight-node  solid  shell  element  in  isoparametric  coordinates:  Sampling  points 
for  ANS  interpolations  for  transverse  shear  strains  (A,  B , C,  D)  and  for  transverse  normal 
strain  ( E , F,  G,  H ). 


To  include  the  constant  stress  in  the  element  (e),  the  orthogonal  condition  of  EAS 
must  hold  in  (3.4),  that  is 

[ S : EdV  = 0 . (3.79) 

We  define  the  following  component  forms  of  the  enhancing  strain  tensor  E as 

E = (£)  0 & (£)  = %iGk  (0)  0 Gl  (0)  , (3.80) 

where  the  enhancing  strain  components  with  respect  to  the  covectors  Gl  (£)  at  any  ar- 
bitrary point  £ are  denoted  by  Eij,  while  those  with  respect  to  the  covectors  Gk  (0)  at 
the  element  center  £ = 0 by  Eki.  From  (3.80),  Ei3  can  be  expressed  in  terms  of  Ekt  as 
follows 


E^  = Eki 


Gl(Z)-Gk( 0)  Gl  (0)  • Gj  (£) 


(3.81) 


where  the  covector  Gk  (0)  can  be  computed  from  the  vector  G;  (0)  by 


Gk( 0)  - GqGi  (0) 
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with 


[Gou]  and  Gold  — Gk  (0)  • Gi  (0)  , 


The  matrix  form  of  (3.81)  is 


(3.82) 


(3.83) 


where  the  components  of  the  enhancing  strain  are  arranged  in  the  same  order  as  in  (3.48), 
and  T0  is  the  matrix  that  transforms  the  strain  components  relative  to  the  basis  { Gt  (0)} 
to  those  components  relative  to  the  basis  {Gj  (£)}. 

Using  the  Column-matrix  form  {5iJ}  for  the  stress  components,  as  in  (3.55),  and 
using  (3.83),  we  can  rewrite  (3.79)  as 

I {Si>}T{ElJ}dV=  J {Sij}Tr0  jlyjdV  = 0 . (3.84) 

B(0e)  B(0e) 


For  constant  stresses,  we  have  the  following  condition  on 


J [Ei^dV  = Jd?ded?  = 0 , (3.85) 

& 

where  J is  the  determinant  of  element  Jacobian  matrix  of  the  mapping  from  the  isopara- 
metric space  □ to  the  initial  configuration  B^  of  element  (e).  Let  j be  defined  by 
using  the  interpolation  matrix  M and  the  element  parameter  a(e)  as  follows 


{£y}  = 2m  ({)«<*>. 

Substituting  (3.86)  into  (3.83),  the  enhancing  strain  j El3  j can  be  written  as 
{^ij}  = Q (£)  a.  , with  Q — ~T0M  . 


(3.86) 


(3.87) 


The  interpolation  matrix  M should  be  constructed  to  satisfy  (3.85)  for  arbitrary 
matrix  The  selection  of  M is  not  unique.  In  the  present  solid-shell  element,  we  use 
sever  parameters  of  EAS  to  pass  the  patch  test  as  shown  below.  To  pass  the  membrane 
patch  test  and  out-of-plane  bending  patch  test,  the  bilinear  polynomials  for  the  transverse 
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normal  strain  E33  are  necessary  (i.e.,  the  minimum  number  of  EAS  parameters  for  E33 
should  be  four.  Therefore,  the  optimal  number  of  EAS  parameters  should  be  seven,  as 
shown  in  the  matrix  M below 


■ f 1 0 0 0 0 0 O' 

0 £2  0 0 0 0 0 

0 0 £2  0 0 0 

. 0 0 0 0 £3  £:£3  £2£3 


(3.88) 


A computationally  more  expensive  choice  for  passing  both  patch  tests  is  to  include  the  tri- 
linear  polynomials  for  E33  and  bilinear  polynomials  for  Eu  (Bischoff  and  Ramm  [1997] 
and  Betsch  and  Stein  [1996]).  In  this  case,  the  number  of  EAS  parameters  is  nine,  with 
the  matrix  M as  shown  below 


' f1  0 0 0 0 0 0 0 O' 

0^2  00000  0 0 

0 0 f1  £2  £J£2  0 0 0 0 

0 0 0 0 0 £3  f1^3  £2£3  £J£2£3 


(3.89) 


The  results  of  our  numerical  experiments  showed  that  there  is  little  advantage  in  using 
(3.89),  since  improvements  compared  to  the  use  of  (3.88)  were  insignificant  (Vu-Quoc 
and  Tan  [2003]).  If  we  enhance  all  the  six  strain  components  [Eu,  2E22,  Eu,  E33,  2E23,  2 E\3], 
the  interpolation  matrix  M contains  complete  sets  of  polynomials  up  to  the  trilinear  one, 
and  thus  corresponds  to  a set  of  30  EAS  parameters  (Andelfinger  and  Ramm  [1993] 
and  Klinkel  and  Wagner  [1997]).  Without  a combination  with  the  ANS  method  to  rem- 
edy the  shear-locking  problem,  the  above  30-parameter  EAS  element  (Klinkel  and  Wag- 
ner [1997])  cannot  pass  the  out-of-plane  bending  patch  test  by  itself  (Vu-Quoc  and  Tan 
[2003]).  Note  also  that  this  set  of  30  parameters  for  enhancing  strain  has  satisfied  the  or- 
thogonality condition  with  the  stress  interpolation  function  which  is  used  in  stress-based 
hybrid  formulation  (Pian  and  Sumihara  [1984]).  Note  also  that  the  new  finite  element  for- 
mulation using  the  above  orthogonality  condition,  but  different  interpolation  functions  of 
stress  and  strain,  has  been  developed  by  Vu-Quoc,  Tan  and  Mok  [2002]  and  Mok,  Tan  and 
Vu-Quoc  [2003],  The  so-called  HybEAS  element  formulation  can  predict  excellent  trans- 
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verse  shear  stresses  through  the  thickness  direction  in  composite  structures  with  elastic 
and  nonlinear  materials. 


CHAPTER  4 

CONSTITUTIVE  MODEL  FOR  SHAPE-MEMORY  ALLOYS 

4.1  Introduction 

Development  of  numerical  tools  to  design  process  of  SMA  devices  for  use  in  three- 
dimensional  structures  has  been  conducted  in  the  last  two  decades.  In  order  to  describe  the 
shape-memory  alloys,  a lot  of  microscopic,  micro-macroscopic,  and  macroscopic  models 
have  been  studied  to  simulate  SMA  devices.  Several  macroscopic  constitutive  models 
that  have  been  able  to  reproduce  the  SMA  material  response  have  been  proposed  in  lit- 
erature, e.g.,  the  exponential  model  by  Tanaka  [1986],  the  cosine  model  by  Liang  and 
Rogers  [1990],  The  modified  cosine  model  by  Brinson  and  Lammering  [1993]  is  one 
of  the  early  numerical  implementations  of  an  SMA  constitutive  model  to  study  the  finite 
element  analysis  of  the  behavior  of  shape-memory  alloys  and  their  applications  using  1- 
D truss  elements.  The  so-called  generalized  plasticity  model  by  Lubliner  and  Auricchio 
[1996]  has  also  been  used  to  simulate  the  shape-memory  materials  response  in  a very  ef- 
fective and  simple  way.  Trochu  and  Qian  [1997]  presented  a phenomenological  material 
law  based  on  plasticity-like  constitutive  model  to  simulate  superelasticity,  one-way  and 
assisted  two-way  shape-memory  effects. 

The  unified  thermodynamical  model  proposed  by  Lagoudas,  Bo  and  Qidwai  [1996] 
has  generalized  to  include  the  exponential  model,  the  cosine  model,  and  polynomial 
model  into  their  thermodynamic  model.  A series  of  papers  (Qidwai  and  Lagoudas  [2000a]; 
Qidwai  and  Lagoudas  [2000/?])  on  the  three-dimensional  numerical  implementation  of 
their  unified  thermodynamic  model  has  showen  that  this  model  can  be  efficiently  used  in 
a wide  area  of  shape-memory  alloys  application  and  has  described  the  different  algorithm 
and  performance  between  implicit  and  explicit  Euler  rule  of  stress  integration  scheme. 
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However,  their  thermodynamical  model  has  some  disadvantage  of  the  model  because  of 
consideration  of  the  change  in  material  properties  and  assumed  only  linear  transformation 
hardening  during  the  phase  transformation. 

The  well-known  return-mapping  algorithm  based  on  implicit  backward  Euler  nu- 
merical integration  scheme  of  the  constitutive  relations  has  been  implemented  by  Au- 
ricchio  [1995],  Govindjee  and  Kasper  [1997],  Terriault,  Trochu  and  Brailovski  [1999], 
Qidwai  and  Lagoudas  [2000c/],  etc.  In  Rengarajan,  Kumar  and  Reddy  [1998],  a return- 
mapping algorithm  based  on  explicit  Euler  integration  scheme  is  developed  based-on  the 
micro-mechanical  SMA  constitutive  model  developed  by  Sun  and  Hwang  [1993a].  More 
general  discussion  on  thermo-mechanics  and  transformation  surfaces  of  polycrystalline 
NiTi  shape  memory  alloy  material  simulated  using  the  unified  thermodynamic  model  are 
presented  in  Qidwai  and  Lagoudas  [2000/?].  Various  types  of  transformation  functions 
such  as  J2,  h — J2,  J2  — or  /1  — J2  — J3  type  are  described  and  compared  with  each 
other  to  explain  the  pressure-dependency  and  asymmetry  under  tension  and  compression 
loading  (Auricchio  [1995];  Gillet,  Patoor  and  Berveiller  [1998]).  More  recently,  Liew, 
Kitipornchai,  Ng  and  Zou  [2002]  has  used  a similar  constitutive  model  to  the  unified  ther- 
modynamical model  with  the  three-point  and  four-point  bending  effects  of  shape-memory 
alloys.  For  more  general  reviews  of  modeling  for  shape-memory  alloys,  readers  refer  to 
Birman  [\991b]  and  van  Humbeeck  [2001], 

The  outline  of  this  chapter  is  given  as  follows.  In  section  2,  Helmholtz  free  energy 
consistent  with  the  principles  of  thermodynamics  is  assumed  to  be  the  sum  of  elastic 
part,  the  chemical,  and  transformation  hardening  part.  A macroscopic  phenomenological 
SMA  constitutive  model  based  on  the  unified  thermodynamical  model,  but  does  not  take 
into  account  the  change  in  material  properties  to  have  the  enhancement  of  numerical 
integration  schemes,  is  also  discussed.  In  section  3,  the  return-mapping  algorithm  to  be 
implemented  into  the  solid-shell  finite  element  formulation  is  also  presented  and  shows 
that  because  of  assumption  of  negligible  changing  in  material  properties,  the  proposed 
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numerical  scheme  shows  a special  case  of  the  closest-point  projection  method.  Finally, 
some  benchmark  cases  such  as  uniaxial  loading  in  symmetry  and  asymmetry  in  tension 
and  compression  as  well  as  three-point  bending  problem  subjected  to  loading/unloading 
in  isothermal  superelasticity  are  presented.  The  proposed  constitutive  model  shows  an 
efficiency  of  numerical  integration  compared  with  the  previous  unified  thermodynamical 
model. 


4.2  Continuous  Time  Constitutive  Model 

SMA  material  behavior  was  modeled  using  the  thermodynamics  of  irreversible  pro- 
cess with  internal  state  variables.  In  this  formalism,  the  constitutive  relation  of  SMA  is 
given  by  its  free  energy  function  and  the  transformation  function  is  obtained  by  intrinsic 
dissipation  rate  during  phase  transformation  flow.  Here  we  use  the  small  strain  additive 
decomposition  hypothesis: 


e = ee  + ept  (4.1) 

where  ee  is  elastic  strain  tensor  and  ept  is  the  phase  transformation  strain  tensor. 

Remark  4. 1 . In  this  thesis,  it  is  assumed  that  at  very  low  deformation  rate  (quasi- 
static process),  the  specimen  behaves  like  an  elastic-plastic  body  in  an  isothermal  process. 
If  the  deformation  rate  is  fast  enough  to  be  truly  adiabatic,  i.e.,  there  is  an  insufficient 
time  for  significant  heat  transfer  from  the  specimen,  the  present  model  will  exhibit  no 
apparent  rate-dependence.  At  increasing  deformation  rate,  the  process  becomes  non- 
isothermal  because  the  heat  produced  by  dissipation  cannot  be  transferred  completely  to 
the  surrounding  system.  At  higher  deformation  rate,  the  process  runs  nearly  adiabatically. 
Rate-dependence  of  Nitinol  taking  into  the  constitutive  model  are  reported  by  Graesser 
and  Cozzarelli  [1994]  and  Juhasz,  Andra  and  Hesebeck  [2000],  Both  models  are  based 
on  visco-plastic  constitutive  model.  | 
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4.2. 1 Free  Energy 

The  Helmholtz  free  energy  is  chosen  herein  for  thermodynamic  potential,  instead 
of  the  Gibbs  energy  for  which  experiments  are  conducted  under  constant  stress  and  tem- 
perature, to  make  it  comparative  with  the  classical  plastic  problem  since  the  Helmholtz 
free  energy  is  the  natural  choice  and  used  widely  when  formulating  constitutive  rela- 
tions of  solids  for  finite  elements  analysis.  Assuming  linear  thermoelastic  response  and 
transformation  hardening  behavior,  the  total  specific  free  energy,  ept , £),  of  a 

polycrystalline  SMA  is  assumed  to  be  equal  to  the  sum  of  the  free  energy  e , T,  e pt ) 
plus  the  free  energy  of  phase  transformation,  the  chemical  free  energy  AjFc/l(T,  £),  and 
the  free  energy  due  to  transformation  hardening  .Fpf(£)  (Huo  and  Muller  [1993],  Brinson 
and  Huang  [1996],  Liew  et  al.  [2002],  Lagoudas  et  al.  [1996],  Rajagopal  and  Srinivasa 
[1999]):2 

H € , T,  ep\  £)  :=  Tel{  e,T,  ept)  + A Tch(T,  £)  + F*(£)  (4.2) 

The  specific  Helmholtz  free  energy  of  thermo-elastic  materials,3  Tel{  e.T,  e pt),  is 
assumed  to  be  (Fung  [1965]) 

£<“)  :=  T(e  - e")  :D:  (e  - e*) 

- i(e  - £*)  : /3(T-T„)  + Q,{ln(Ar-[r-r„]}(4.3) 

2 For  example,  in  one-dimensional  plasticity  problem  we  can  write  the  strain  energy  function  W as 
follows: 

T(e,  ep,  a)  = ep)  + a ) = ~(e  - tp)E{t  - ep)  + | a Hlso  a 

where  HlS0(  a ) is  nonlinear  isotropic  hardening  function  and  a is  an  internal  variable  with  which  the  plastic 
strain  can  usually  be  replaced  for  one-dimensional  analysis.  From  which  we  could  derive  the  constitutive 
relations  by  derivative  with  respective  to  each  of  the  arguments. 

3 In  thermodynamics,  intensive  properties  which  is  defined  per  unit  mass  (also  called  specific)  are  used. 
On  the  other  hand,  in  solid  and  structure  mechanics,  the  quantities  per  unit  volume  are  widely  used,  e.g., 

£ = G + pTrj  + a : e 

where  £ , Q,  and  rj  are  internal,  Gibbs  energy  and  entropy  per  unit  volume.  Here  we  follow  the  conventions 
used  in  the  thermodynamic  principle  since  SMA  behavior  is  analyzed  by  thermo-elasticity. 
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where  (3  is  referred  to  the  thermal  stress  coefficient  tensor,  defined  as  (3  = D : a where 
a is  the  thermal  expansion  coefficient  tensor.  Sometimes  the  elastic  strain  energy  Tel  is 
called  the  mechanical  energy  Fnech . c0  is  the  specific  heat  at  reference  state. 

The  change  in  chemical  free  energy  due  to  phase  transformation  is  a material  prop- 
erty which  is  a function  of  temperature  and  entropy.  Different  phases  of  the  same  mate- 
rial have  different  chemical  free  energies  (Huo  and  Muller  [1993],  Lu  and  Weng  [1997], 
Huang  and  Brinson  [1998],  Auricchio  and  Sacco  [2001],  Qidwai  and  Lagoudas  [2000a, 
eq.21],  Raniecki  and  Lexcellent  [1994]). 


()  .=  rH-  Tt  :=  £„  - Tm 


(4.4) 


where  rj0  is  the  entropy  defined  as  follows: 

Vo  ■=  Vo  + Z(Vo  ~ Vo)  = Vo+  £At7o  (4.5) 

£o  :=  £o  + £(C  ^o  ) • = ^o  "T  £A£q 


where  Ar]0  :=  rj ft1  — t]q  and  A£0  £(fr  — £$.  Chemical  forces  arise  from  the  difference 
of  free  energy  between  austenite  martensite  and  act  as  a driving  force  of  the  transforma- 
tion. In  order  to  induce  the  austenite  phase  to  martensite  (A  — > M)  transformation,  the 
chemical  free  energy  of  the  martensite  phase  must  be  lower  than  that  of  austenite  phase 


(Meyers  and  Chawla  [1999,  p.498];  Ortin  and  Planes  [1988];  Gil  and  Planell  [1999]). 

The  specific  free  energy  for  kinetics  and  hardening  during  tranasformation  is 

assumed  to  be  a quadratic  function  in  f (Huo  and  Muller  [1993];  Lagoudas  et  al.  [1996]).4 
For  example,  for  the  isotropic  transformation  hardening  case. 


•F*(0  :=  =rK £2  + K 


'?  + iexp( -/?"?)) 


(4.6) 


where  k is  constant  material  parameter  to  be  determined  or  given,  representing  the 
isotropic  hardening  modulus  and  u = A,  M for  austenite  and  martensite  phase,  respec- 
tively (For  more  detailed  description,  see  section  4.2.6). 

4 Huo  and  Muller  [1993]  proposed  the  mixing  free  energy  of  phases  due  to  their  coherence:  Tmix  = 
M - A£?  where  A is  a non-negative  constant,  called  the  coherence  coefficient. 
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We  now  can  obtain  the  constitutive  relations  from  the  above  thermodynamic  poten- 
tial as  follows: 


dT 

Pj^  = D : e - D : ept  - (3  [T  - T0]  = a 


(4.7) 


P 


dT 
deP l 


= -D  : e + D : ept  + /3[T  -T0]  = - cr 


From  (4.7)  and  (4.8),  we  can  show  that 


dT 


dT 


Pde  P d ept 


(4.8) 


(4.9) 


which  is  very  similar  to  the  case  of  plasticity. 

4.2.2  Internal  Entropy  Production  Rate 

Based  on  the  principles  of  thermodynamics,  we  know  that  the  specific  Gibbs  poten- 
tial is  related  to  the  specific  internal  energy,  S(cr  , rf)\ 


£(a,rj):=Q  + Tr]+-a  : e 

P 

where  e is  the  strain  tensor.  The  Gibbs  potential  is  thus  given  by 

Q :=  T cr  : e = £ — r/T cr  : e 

P P 

where  T is  the  specific  Helmholtz  free  energy,  defined  by 


(4.10) 


(4.11) 


T:=£-pT 


(4.12) 


from  which  the  expression  of  differentiation  in  time  for  £ is  given  by 


• • 


£ = T + r]T  + VT 


(4.13) 


From  the  first  law  of  thermodynamics 


• • 

p£  = cr  : e + Q = cr  : e+pr  — div(q) 


(4.14) 
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where  r the  internal  heat  source  per  unit  mass  and  q is  the  heat  flux  vector  per  unit  area 
through  the  surface  of  the  body.  The  second  law  of  thermodynamics  is 


PV  - —pr  + div 


> 0 


(4.15) 


Combining  both  the  first  and  the  second  law,  the  Clausius-Duhem  inequality  can  be  writ- 
ten as: 


p[Tfi  - £}  + or  : e - ^grad(T)*q  > 0 


(4.16) 


Using  (4.13)  and  (4.16) 


• • 1 

— pT  — pr/T  + cr  : e — — grad(T)*q  > 0 


(4.17) 


Since  T = T(  e , T,  £,  ept)  and  the  total  differentiation  in  time  gives 


Tie  TP  ept)  = — • e + — T + — £ + ^ ■ e pt 
1 J de'e  + dT1  + dC  + de*  ' 


(4.18) 


df 

• 

dT 

a ~ 

o e 

: e - p 

V+  dT. 

T 


~PJ7*  : lPt  ~ ~ rad(T)-q  > 0 


(4.19) 


pt  7, 

The  total  entropy  production  rate  can  be  decomposed  into  two  parts: 

pT  % = PTim  + PT  7 h (4-20) 

Here  we  assume  the  thermal  term  is  to  be  always  positive  (Malvern  [1969,  p.268]) 

PTlh  '■=  -^gfad(r)*q  > 0 (4.21) 


which  is  called  heat  or  thermal  entropy  production  rate,  which  means  that  the  direction 
of  grad(T)  is  opposite  of  q,  i.e.,  the  heat  flows  from  hot  to  cold  region.  Mathematically, 
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the  angle  between  q and  grad(T)  is  greater  than  right  angle.  We  arrive  at  the  Clausius- 
Duhem  inequality: 


dT' 

• 

dT' 

a ~P^~ 
de 

: e -p 

V+dT 

T-p 


dT 
d ept 


v 


pT'lm 

and  finally  the  rate  of  local  intrinsic  dissipation  can  be  defined  as  follows: 


(4.22) 


pTirn--=-Pj^t--^-Pi^>  0 (4'23) 

Using  the  constitutiion  relation,  the  internal  entropy  production  rate  can  be  expressed 

PTim  ■=  d : evt  - p— I > 0 (4.24) 

Note  that  for  pure  elastic  processes  (4.9)  is  necessary  and  sufficient  for  (4.22),  while  for 
elastic-transformation  processes  (4.9)  together  with  (4.24)  are  only  sufficient  for  (4.22). 
4.2.3  Constitutive  Relations 

Since  dissipation  is  associated  to  phase  transformation,  the  entropy  production  must 
vanish  with  £ so  that  (4.22)  give  the  elastic  laws  and  entropy  production  inequality.  Equal- 
ity in  4.22  holds  for  all  values  of  e and  T.  That  is,  the  following  set  of  constitutive 
relations  can  be  obtained  by  identically  satisfying  the  inequality 

dT  J dT 

a :=  pd7  and  77  :=  ~df  (4,25) 

Using  the  expression  of  Helmholtz  free  energy  previously  defined  (4.2),  we  obtain  the 
following  relations: 


a :=  D : 

e - p [T  - T0]  - D : ept 

= D : 

(e  -D-1  : /3[T  — T0]-  e^) 

= D : 

(e  - a [T  — T0]  — ept) 

= D : 

(e  - e9  - ept ) 

(4.26) 
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where  e 0 is  the  thermal  strain  and 


V := 


T 

(3  + c0ln(— -)  + r]0 
Jo 


(4.27) 


Remark  4,2.  For  the  case  of  isotropic  thermo-elastic  materials,  the  stress-strain  re- 
lationship is 


a — ^kkfiij  + 2/uejj  — (3A  + 2 n)a8ij[T  — To] 

— + 2 — 3 A aSij[T  — To] 

= XtkkSij  + 2fieij  - 05^  - T0]  (4.28) 

where  the  bulk  modulus  is  K = — crfcfc/efcfc,  A and  p are  Lame’s  constants.  I 


4.2.4  Transformation  Function 

The  local  internal  entropy  production  rate  (4.24)  has  been  written  by 

pTym:=  cr  : ept  - p~l  > 0 (4.29) 

The  above  expression  can  also  be  simplified  and  generalized  using  dissipative  power  as 
follows  (Malvern  [1969,  p.268];  Maugin  [1992,  p.44];  Tanaka  and  Iwasaki  [1985]): 

0(X)  :=  Y : X > 0 (4.30) 


where 


and  (4.3  n 

where  Y is  called  the  generalized  irreversible  thermodynamic  force  and  X is  the  general- 
ized thermodynamic  flux.  Using  the  rule  of  transformation  strain  flow,  ept  = which 
will  be  discussed  later,  the  local  internal  entropy  production  rate  can  also  be  rewritten 


pT'frn  '■  = 


5$ 

1 

q 

Qj 

~ P^7 

£ = n > o 


(4.32) 
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Inequality  (4.32)  is  the  thermodynamic  constraint  to  the  SMA  material  system.  Here  the 
thermodynamic  force  conjugate  to  £,  Y is  defined  as 


Y :=  a 


d<b  dT 
: da 


(4.33) 


The  transformation  function,  like  the  yield  function  in  plasticity,  is  thus  defined  as 


+Y-Yf,  (|>0) 

-y-y;,  (|<o) 


(4.34) 


where  Y* , will  be  described  and  determined  analytically  later,  is  the  measure  of  in- 
ternal dissipation  during  the  phase  transformation.  During  phase  transformation,  the 
driving  force  to  martensite  transformation  is  countered  by  the  resistance  force.  The  en- 
ergy dissipation  associated  with  the  austenite-martensite  interfacial  friction  and  the  self- 
accommodating  twinning  motion  will  result  in  stress  hysteresis  in  the  superelastic  cycle 
and  interfacial  friction  between  phases  (Sun  and  Hwang  [1993a],  Lagoudas  et  al.  [1996], 
Lu  and  Weng  [1997],  Rajagopal  and  Srinivasa  [1999]).  That  is,  whenever  £ f 0 the 
driving  force  Y must  be  equal  to  the  dissipative  resistance  Y*. 

Thus,  thermodynamic  function  during  forward  phase  transformation  can  thus  be 
rewritten  as 


0 :=  a : — -(-  (pAr,0)  T - (pA£0)  - _ y*  (4.35) 

da 

We  now  decompose  additively  the  transformation  function  into  stress  term  and  tempera- 
ture and  martensite  volume  fraction  term  as  follows:5 


<$>f{a,T,0  :=  [R°f(T,0  + Rhf(T,^]  <0 

= ~ Rf(T,£)  < 0 


(4.36) 


5 In  plasticity  of  1-D  model,  the  yield  function  can  be  written  as 

fy(a,  a ) = |er|  - [cr0y  + Hlso  a ] < 0 

where  Hiso  is  the  linear  isotropic  hardening  modulus  and  a represents  the  evolution  of  isotropic  hardening. 
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where 


$(0-) 


(4.37) 


Rf(T,t) 


R°f(T^)  + Rhf(T^) 


(4.38) 


which  represents  the  radius  of  circle  of  von-Mises  type  transformation  function. 


Remark  4.3.  The  present  model  can  easily  be  extended  to  cover  the  plastic  zone, 
beyond  the  fully  detwinned  region  as  follows: 


$/(cr) 


fl/CT.f)  + Rl(T,0\ 


Gy  + k(ot) 


< 0 


(4.39) 


where  a is  isotropic  hardening  evolution  parameters  and  <7°  is  the  initial  plastic  yield 
stress.  Following  the  classical  plasticity,  plastic  flow  rule  can  be 


<9$ 

• j)  • 

6 =7^ 


where  7 is  the  plastic  multiplier  parameter  (Ford  and  White  [1996]). 


(4.40) 

I 


The  transformation  function  of  stress-related  is 

« <94>  / — 

®f{v)  :=  ^ = V3J2eL  (4.41) 

Note  that  the  above  function  is  a function  of  only  in  stress  to  have  the  associative  transfor- 
mation flow  rule  and  later  to  preserve  the  closest-point  projection  schemes  in  time  discrete 
model. 


R°JT,  0 :=  - (pAVo)  T - (pAS0)  - Y*f 


which  plays  a role  of  initial  yield  stress  in  plasticity  and 


Rhf(,T,t)  :=  + 


d£ 


(4.42) 


(4.43) 
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which  is  the  evolution  of  the  isotropic  hardening  in  plasticity. 

The  thermodynamic  function  during  reverse  phase  transformation  can  thus  be  rewrit- 
ten as 

$r(cr,T,0  :=  -a  : — - (pAri0)T  + (pA£0)  + -^ Y*  (4.44) 

Decomposing  additively  the  transformation  function  into  stress  term  and  temperature  and 
martensite  volume  fraction  term  yields  as  follows: 


:=  *r(<r)-  R°r(.T,t)  + Rhr(T,Z) 


(4.45) 


where 


flr(T.f) 


'A?(r,0  + #(T,f)] 


(4.46) 


which  represents  the  radius  of  circle  of  von-Mises  type  transformation  function.  The 
transformation  function  of  stress-related  is 


cr ) :=  - a : — = -yJzI2eL 


(4.47) 


R°r  (T,  0 :=  + [(pA?7o)  T - (pA£0)  + Y*]  (4.48) 

playing  a role  of  initial  yield  stress  in  plasticity  and 

*?CT,0:=-- IP  (4,49) 

which  is  the  evolution  of  the  isotropic  hardening  in  plasticity. 

The  above  additive  decomposition  of  transformation  function  has  great  flexibility 
of  implementation  by  replacing  it  with  other  types  of  transformation  functions  such  as 
^((j)  = $(J2),  4>(cr)  = <!>(/!,  J2),  or  4>(cr)  = <!>(J2,  J3)  with  which  the  classical 
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plasticity  has  an  analogy.  In  this  section,  we  review  briefly  a few  different  types  of  trans- 
formation functions.  Assume  that  the  general  transformation  function  can  be6 

bJ3  V 


$(<r)  :=  eLJ3J2  1 + 


(4.50) 


W2) 

Case  1:  When  6 = 0 and  C = 1,  <!(  cr ) = l>(J2). 

When  6 = 0 and  ( = 1,  the  function  can  be  reduced  to  the  well-known  von-Mises  type 
function  of  J2  flow  such  that 


<!>(  a ) <f>(  J2)  (-l\]3J2 


(4.51) 


The  flow  rule  can  be  described 


[pt  -d^a,T,Q 

da  da 


(4.52) 


which  leads  to  the  fact  that  the  direction  of  phase  transformation  strain  tensor  is  the  normal 
to  the  transformation  surface,  proved  as  follows: 

d${a)  d$(J2)  dsj 3T2 


da 


= eL- 


d a da 


<9. 


3 a' 
2~d~ 


(4.53) 


where  the  equivalent  stress  cr  is  defined  as 


1/2 

a :=  •?  -s  : s > = 


cr'||  = \/3  J2 


(4.54) 


Thus,  the  transformation  direction  tensor  A is  obtained  to  have  the  following  form 


d<f> 


3 1 


da  2 o 


(4.55) 


6 The  general  transformation  function  is  based  on  Qidwai  and  Lagoudas  [2000/;],  but  has  a little  dif- 
ference in  that  \/3  in  the  second  term  here  is  dropped  to  be  reduced  to  Gillet  et  al.  [1998]. 
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Case  2:  When  b ^ 0 and  £ = 1/2,  4>(  cr  ) = $(  J2,  J3). 

We  can  also  recover  the  same  form  of  Gillet  et  al.  [1998]  when  C = 1/2  as  follows:7 

/ h 7 \ 

*(<r):=£iV^l  + ^j  (4.56) 

, — / h J \ 1/2 

4>(J2,J3,T,e'»)  :=  A(l  + - *(T,  £*)  = 0 (4.57) 

By  squaring  on  both  sides  of  transformation  function,  we  can  recover  Gillet’s  model  as 
follows: 


4>(<T)  = <f>(j2,j3,r,  e^) 

where 


h 


k2(T,  ept)  = 0 (4.58) 


b 


V"27  a2c  - of 

2 + at 


(4.59) 


where  ac  is  the  initial  yield  stress  in  compression  and  at  is  the  initial  yield  stress  in 
tension. 

Case  3:  When  b ^ 0 and  £ = 1,  $(  cr ) = <!>(  J2,  J3). 

In  (4.50),  letting  £ = 1 for  simplicity  will  give  the  following  formula8 


Her)  := 


J 3 


— 3 ( y + 6-jr  ) — Cl y 3t/2  + e/,\/36— - 

J2/  J2 


(4.61) 


7 In  the  original  Gillet’s  model,  \/3  is  transfered  to  the  left  hand  side,  similar  to  the  classical  plasticity 
where  the  yield  function  is  expressed  either  /y(  a , a")  :=  x/3^2  - o^a")  = 0 or  /y(  <r  , a")  :=  vGJ  - 
= 0. 

x In  elasto-plastic  problem,  similar  type  of  yield  function  of  $(  J2,  J3)  has  been  used  by  Auricchio  and 
Taylor  [1999]  which  is  referred  to  as  Prager-Lode  yield  function: 

fy(cr)  = y/2J^  + P^-  -Gy  < 0 

j 2 


where 


and  = 


2 (Tc  4~  (Tj  3 Gc  + Gt 

where  crc  is  the  initial  yield  stress  in  compression  and  <rt  is  the  initial  yield  stress  in  tension. 


(4.60) 
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Next,  we  want  to  determine  the  derivative  of  <I>  with  respect  to  stress.  Using  the  transfor- 
mation function  above 

d$(<r)  = af (j2,  j3)  _ a$  dj2  a$  dJ3 

da  da  dJ2  da  + dJ3  d a 

=:  Fi(J2,  J3)  + F2(J2,  J3)  (4.62) 


where 


Fi(J2,J3)  := 


d<l  dJ2 


dJ2  da  ’ 

Recalling  some  of  tensor  algebra  in  the  following 

dJ2  d (1 


F2(J2,  J3)  := 


a$  dJ3 
dJ3  da 


'v 

where  the  stress  deviator  is  defined  as 


dati  d*u  \2SijSij  1 ~ Sij 


(4.63) 


(4.64) 


Sij  (Ttj  g ^kk^ij 


(4.65) 


dsn  dan  1 dakk  1 1 

-Stj  = SipSjq  - -SpkSgkSij  = SipSjg  ~ ^S^pg 


dapq  dopg  3 do 


pq 


or 


(4.66) 


as  1 

= 14- -101 
da  3 


For  the  third  invariance  of  stress  tensor. 


1 


J'i  det(  a ) 

the  first  derivative  of  J3  with  respect  to  stress  is  given  by 

aj3 


dot  tlj 


'o 


21 


SikSkj  g 2 SrsSrs°ij  — $ikSkj  2 2^i 


(4.67) 


(4.68) 


(4.69) 
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The  first  derivative  of  <&(  J2,  J3)  with  respect  to  the  stress  is 

a$(J2,  J3)  _ crt^dJ^  <9J3 
dcjij  dJ2  duij  d,J:i  dcrij 


4.2.5  Transformation  Strain  and  Transformation  Flow  Rule 

Here  infinitesimal  strain  hypothesis  is  used.  Small  strain  additive  hypothesis  is  a 
reasonable  assumption  for  polycrystalline  SMA  materials  undergoing  constraint  defor- 
mations (Rengarajan  et  al.  [1998,  p.  1497]).  Thus  apparently  plastic,  large  deformation 
is  mostly  recoverable,  and  hence,  a small  strain  assumption  is  reasonable  as  long  as  the 
application  does  not  involve  finite  rotations. 

e = ee  + ept  (4.71) 

where  e e is  the  elastic  strain  tensor  and  e pt  is  the  phase  transformation.  Since  the  phase 
transformation  is  almost  volume  conserving,  the  following  assumption  can  be  introduced 
to  relate  the  evolution  of  the  phase  transformation  strain,  e pt,  to  the  evolution  of  the 
martensite  volume  fraction: 

= (4.72) 

O <J 

where  £ is  the  volumetric  fraction  of  martensite  phase  and  4>  is  the  transformation  func- 
tion. 

In  case  of  J2  flow,  the  gradient  of  transformation  surface  is  assumed  to  have  the 


following  form: 
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(4.73) 


where  eL  :=  e^ax  is  the  maximum  uniaxial  transformation  strain.  We  denote  the  total 
inelastic  phase  transformation  strain  at  the  end  of  the  A — > M transformation,  i.e.,  trans- 
formation strain  at  the  end  of  forward  phase  transformation,  parametrically  as  which 
represents  the  transformation  strain  due  to  pre-strain  (Lagoudas  and  Bhattacharya  [1998, 
p.335],  Sun  and  Hwang  [1993a],  Rengarajan  et  al.  [1998]).  Note  that  the  value  of  this  ten- 
sor will  depend  on  the  extent  to  which  the  partial  or  full  A — *•  M transformation  has  taken 
place  in  the  previous  transformation  cycle.  As  a result,  during  the  reverse  transformation, 
the  stress  applied  to  the  opposite  direction  of  the  accumulated  transformation  strain  for 
either  fully  or  partially  transformed  martensite  will  help  the  reverse  transformation.  The 
equivalent  transformation  stress  is  defined  as  (Malvern  [1969,  p.364]) 


where  J2(  <x ) is  the  second  invariant  of  stress  deviator 


s = <*  ~ ^(trtr)  1 


(4.75) 


The  equivalent  transformation  strain  is  defined  as 


(4.76) 
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where  ept  is  the  transformation  strain  at  the  reversal  of  phase  transformation  (A  —>  M) 
and  J2(  ept)  is  the  second  invariant  of  transformation  strain  tensor.9 

Remark  4.4.  The  form  of  e pt  is  satisfactory  for  the  case  of  proportional  load- 
ings, but  not  for  taking  account  of  the  non-proportional  loadings  for  reorientation  effect. 
Macroscopically  proportional  loading  of  shape-memory  materials  generally  results  in  mi- 
croscopically non-proportional  loading  of  the  martensite  phase  in  SMAs  and  austenite  in 
SMAs.  During  non-proportional  loading,  reorientation  may  accompany  transformation 
(Boyd  and  Lagoudas  [1996,  p.807],  Auricchio,  Taylor  and  Lubliner  [1997]).  Juhasz  et 
al.  [2000]  have  used  the  model  for  the  case  of  non-proportional  loading  and  isothermal 
pseudo-elasticity  to  include  reorientation  effects. 

ept=A£+A£  (4.77) 

where  the  first  term  on  the  right-handed  side  means  the  change  of  ept  with  the  phase 
transformation,  the  second  term  means  the  change  of  the  stress  direction  due  to  inhomo- 
geneous stress  fields  where  induces  a reorientation  effects  in  SMAs.  I 


The  normality  of  the  total  transformation  strain  rate  can  be  found10 

-L  (l>o) 


d<bf . 


ept  = 


(4.78) 


where  the  gradient  of  transformation  function  of  von-Mises  type  J2  flow,  $(<x,T, £) 
with  respect  to  stress  tensor  yields 


3<f> 

da 


+ A , (|  > 0) 

-A,  (£<  0) 


(4.79) 


which  shows  the  normality  to  the  transformation  function  ($)  in  stress  space  and  the  flow 
rule  for  the  transformation  strain  evolution  to  be  associated  in  the  essence  of  J2-plasticity. 

9 See  Appendix  for  derivation  of  1-D  case  recovery  from  this  3-D  formula. 

10  See  Appendix  for  detailed  description  of  principle  of  maximum  transformation  dissipation  for  nor- 
mality and  associative  flow  rule  in  stress  space. 
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•Pj 


Figure  4. 1 . Associative  transformation  flow  rule  and  its  normality. 

Remark  4.5.  It  is  noticed  that  £ acts  like  a transformation  parameter  in  analogy  to 
the  plastic  parameter  7 in  associative  plasticity. 


ep  = 7 dfy 


dcr 


(4.80) 


where  ep  is  the  plastic  strain  rate  and  fy  is  the  yield  function  defined  as  fy(cr . otv)  = 
J2  — k(  a ")  for  the  von-Mises  yield  criterion  where  a " is  the  internal  variable  for  the 
strain  hardening  and  J2  is  the  second  invariant  of  stress  deviator.  I 

4.2.6  Transformation  Hardening  Function 

The  free  energy  Jrpt(£)  = //l(£),  is  responsible  for  the  kinetics  of  transformation 
flow  and  the  transformation-induced  strain  hardening  during  phase  transformation.  The 
equations  of  line  of  stress  in  Figure  4.2  read  as  follows: 


pM  — -M 
r2  ■—  °2 


CM(T  - Mf)  = 0 and  F?1  :=  af  - CM (T  - Ms)  = 0 


F2A  :=a£ -CA(T-Af)  = 0 and  FA  :=  aA  - CA(T  - As)  = 0 (4.81) 

so  that  at  T = T (here  T = Af),  a?1  = a™,  a\l  = a*1,  af  = aA,  and  a£  = a£.  By 
analogy  to  the  well-known  Voce’s  exponential  strain  hardening  (Voce  [1948]),  we  have 


M 


_M 

°0  ■—  eL 


-M  -M 
aoo  ~a0 
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6LCm  (Ma  - Mf)  = -(pAp0)  (Ms  - Mf ) 


(4.82) 


aoo  a0 


tL 


-A 

aoo  a0 


— £lCa  (Af  — As)  — —(pAp0)  (Af  — ^4S) 


(4.83) 


where  a ^ is  called  the  saturation  stress  and  a^f  is  the  initial  stress  at  the  onset  of  phase 
transformation.11  Finally,  we  can  assume  the  following  transformation  hardening  func- 
tion which  includes  both  linear  hardening  and  Voce-type  exponential  nonlinear  hardening 
function.  When  / 3M  = (3A  = 0,  the  function  can  be  reduced  to  linear  function  which  has 
been  extensively  used  in  Lagoudas  et  al.  [1996]. 


1 


fn(0  ■= 


i /—AT  \ 

+ (aoo  ~~  a0  ) 


M \ 


£ + ^exp(-/?"0 


1 


+mi£ 

1 


(pAp0)(Ms  - Mf)? 


- (pAVo)  (Ms  - Mf) 

+mi£ 


7+4_exp(-/3"0)-4. 


(4.84) 


/*( f)  := 


- *S) 


1 


« + ^exp(-n))-V 


:=  -^(p^Vo)(Af  - As)£2 


- (pAp0)  (Af  - As) 


f + ijexpt-/^)  )-± 


See  section  4.4  for  detailed  description  about  Clausius-Clapeyron  equation  and  pArjo  = —CMei. 
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Figure  4.2.  Stress-temperature  (a  — T)  phase  diagram  of  superelasticity. 


-mi£ 


(4.85) 


where  mi  is  used  for  continuity  at  £ = 1,  i.e.,  fff{  1)  = fA (1). 


a/f(0 


- ^o')  £ + (*«  - -oM)  (l  - exp(-/3A/0)  +mi 

S" v* ''  S — v*  y 

linear  hardening  nonlinear  hardening 

- (pArio)  ( Ms  - Mf)£  - (pAp0)  (Ms  - Mf)  (l  - exp (-/3A'0) 


+TOl 


(4.86) 


df£(  0 


(^>  ~ ) (!  - exp(-/?^0)  - 

- (pAp0)  (Af  - As)£  - (pApo)  ( Af  - As)  (l  - exp(-/^0) 

-mi  (4.87) 


Nonlinear  isotropic  hardening  is  characterized  by  the  hardening  exponent  /3M  and  f3A.  If 
pM  _ pA  _ tjien  we  recover  iinear  hardening  type.  The  following  properties  are 
assumed  to  be  preserved  for  /),(£): 
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• The  parent  austenite  phase  is  stress-free  if  no  external  mechanical  loading  is  ap- 
plied. This  condition  is  satisfied  by  selection  fh{£)  to  be  zero  at  fully  austenite 
state,  i.e.,  fh( 0)  = 0,  i.e.,  f£( 0)  = fj*( 0). 

• Since  the  partial  derivative  of  //*(£)  in  transformation  function  as  the  onset  of  phase 

dfM(  0)  dfA(0) 

transformation,  --4^  — = ^ ■ = 0.  That  is,  in  elastic  range,  there  is  no  effect  of 

phase  transformation. 

• The  function  fh{0  must  be  continuous  during  the  phase  transformation  as  well 
as  at  return  points  where  the  phase  transformations  direction  is  characterized  by  a 
change  in  the  sign  of  £.  Here  mi  is  used  to  enforce  the  continuity  at  £ = 1,  i.e., 
fhi1)  = fh(l)  (Lagoudas  et  al.  [1996]). 

• The  function  //>(£)  must  be  non-negative.  This  requirement  can  be  satisfied  by 
appropriately  choosing  material  constants,  i.e.,  /^(£)  > 0 and  f£{£)  > 0. 

• The  range  of  £ is  given  by  0 < £ < 1 (Figure  2.1). 


Applying  the  above  constraints,  a form  for  the  function  /*(£)  can  be  selected  as  follows: 


A(0 


/no,  (£>o) 
/.A(0,  (£  < o) 


(4.88) 


See  Figures  4.3  and  4.4  for  the  profiles  of  the  assumed  hardening  function  which  are 
satisfied  above  properties.  Figures  in  the  first  row  are  for  the  linear  hardening  function, 
figures  in  the  second  row  are  for  the  nonlinear  hardening  function,  and  figures  in  the  third 
row  are  for  the  superposition  of  linear  and  nonlinear  transformation  function. 

4.2.7  Transformation  Consistency  Condition 

In  order  to  obtain  the  evolution  equation  for  the  internal  state  variable  f , we  use  the 
consistency  condition.  The  constraints  on  the  evolution  of  the  martensite  volume  fraction 
are  expressed  in  terms  of  the  Kuhn-Tucker  conditions  (Simo  and  Hughes  [1998,  p.77]) 
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Transformation  hardening  function  and  its  derivative 


Figure  4.3.  Transformation  hardening  function  fh(£)  and  its  derivative  when  (3  = 
0.01  ss  0. 


Transformation  hardening  function  and  its  derivative 


Figure  4.4.  Transformation  hardening  function  fh(£)  and  its  derivative  when  (3  — 

1.0. 
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as 


12 


Forward,  (A  — > M)  =>  | > 0,  $(  a , T,  £)  < 0,  £<J>  = 0 
Reverse,  (M  — > A)  ==$■  £ < 0,  <3>(  cr  , T,  f ) < 0,  £<3>  = 0 


(4.89) 


Here,  $(er,T,  £)  signifies  the  transformation  function  that  defines  the  elastic  domain 
such  as  in  yield  surface  in  plasticity.  The  inequality  constraints  on  <h(  a , T,  £)  may  also 
be  called  as  the  transformation  consistency  condition  and  described  as  follows: 

1.  Elastic  response:  For  <f>  < 0 and  t>  < 0 ==>  £ = 0 

2.  Forward  transformation  {A  — > M):  £ > 0 =>•  $ = 0 and  $ = 0 

• • 

3.  Reverse  transformation  (M  — > A):  £ < 0 =>•  $ = 0 and  <f>  = 0 

With  noting  the  fact  that  there  is  no  change  in  transformation  strain  during  transfor- 
mation, we  have  the  following  consistency  conditions. 


1 . During  the  forward  transformation  (A  — > M), 


• <9$  f <9<f>  f • • 


dT 


The  evolution  equation  for  is  thus  given  by 


<9$/  ^ + 


tf  = 


da 


dT 


3<f> 


with 


/ 


dtf 


d<S>f 

^7 


^0 


2.  During  the  reverse  transformation  (M  — > H), 


$r(cr,T,£r)  = 


d$r 

da 


<9<f>r  • <9$,  • 


a 


— -T  + — -e  = 0 
dT  d£r 


The  evolution  equation  for  fr  is  thus  given  by 


d$r 


■ a + 


^Lf 

dT 


d$r 

d£r 


d$r  , „ 

with  — - — ^ 0 

dir 


(4.90) 


(4.91) 


(4.92) 


(4.93) 


12  The  constraints  are  derived  from  the  principle  of  maximum  transformation  dissipation  which  is  sim- 
ilar to  the  principle  of  maximum  plastic  dissipation.  See  Appendix  for  more  detailed  description. 
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4.2.8  Continuum  Tangent  Moduli 

In  an  incremental  displacement-based  finite  element  analysis,  the  increment  of  the 
stress  obtained  through  the  implementation  of  the  incremental  SMA  constitutive  model 
for  given  increments  of  strain  and  temperature.  The  SMA  constitutive  model  can  be 
written  in  an  incremental  form  as  follows: 


dcr  — L\  de- 1-0  dT 


(4.94) 


where  C is  the  continuum  tangent  material  moduli  tensor  (CMAM)  and  0 is  the  con- 
tinuum tangent  thermal  moduli  tensor  (CTEM).  These  two  moduli  are  needed  for  the 
execution  of  the  global  finite  element  analysis. 

In  order  to  derive  C and  0 , the  constitutive  relation  is  recalled  (4.26) 


a :=  p—  = D : (e  - a [T  — T0]  — e*)  (4.95) 

Differentiating  with  respect  to  time  (4.95),  writing  it  in  incremental  form  and  noting 

(4.96) 


<9$ 

dept  = 


dcr 


the  incremental  form  of  the  stress  yields 


dcr  = Id  : < de  — a dT  — 


Ja 


dd 


(4.97) 


During  Forward  Transformation  (£  > 0). 

First,  for  the  forward  transformation  case,  using  (4.79),  the  above  equation  reduces 
to 


dcr  = D : I de  — a dT  — — da 

dcr 


(4.98) 


Recalling  the  consistency  condition  (4.90)  = 0 is  obtained  by  differentiation  of  the 

transformation  function  such  that 


<9$  <9<f> 

d$(cr,T,a)  = ^r:dcT+—dT+—da  = 0 


(4.99) 
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By  substituting  d a from  (4.98)  into  (4.99),  i.e., 


<9$  f «9<f>  ] <9$ 

_:D  :)de.adT._dA  + _dT+_d(=0 


(4.100) 


da 


: D : d e 


<9$ 


<9T  (9  a 


: D : a.  } dT 


d$  _ d<f>  <9$ ) , 

d^:D'd^~d(\di=0 


(4.101) 


Arranging  for  d£  results  in 


dt  = 


da 


^ _ l ' 

:D:de+  sr-^;D:QrT 


where 


f _ <9$  d<f> 

a :=  < - — : D : 


(4.102) 


(4.103) 


da  da  <9£  j 
Now  (4.102)  can  be  used  to  eliminate  d£  in  (4.98)  and  obtain  the  relationship  between 


stress,  strain  and  temperature  increments  as  follows:13 

<9$ 


da  = D < d e — ctdT 


da 


(4.104) 


= D : d e — D : a dT  —ID: 


= D : d e — D : a dT 


da 


dd 


»S)(S— )H 

D,iy 


= D : de 


o cr  y \<9T  da 

)(»»")(£— 


-dT 

i 

1 

a 


(4.105) 


-D  : a dT 


13 


See  Appendix  for  detailed  tensor  algebra  for  the  derivation. 
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D 


+ D 


d<&  \ d$  1 _ 
— -dT 


da  J dT  a 
d$\  ( d§ 


= D-  D 


da  J \da 
d<$>  \ 


D : a)  -dT 
/ a 


(4.106) 


D-  D 


da  J 
<9$ 


D 


\ i 


D 


da 

d$\  d$  1 


D 


d a J a , 
\ 1 


- } : de 


da 


a.  dT 


= D-  D: 


da)  dT  a 
<9$  \ 


dT 


-r 

£ : de 


D : 


da  J 
d$ 


D : 


\ 1 


da  / a 


d e 


da 

+ edT 


' \ 1 
D : - — - : a 
da  / a 


(4.107) 


The  continuum  tangent  material  modulus  tensor,  Cf,  and  continuum  tangent  thermal 
modulus  tensor,  0 /,  for  forward  transformation  are  then  given  by 


£f  = D — A®  A 


(4.108) 


where 


©/  = -C,  : a-~~A 
dT  y/a 


\ d$>  d$  34> 
a < — — - : D : 


da 


, / d$\  1 

o Pc  ( and  A :=  D : — —= 

da  dl;  ) \ da  J y/a 


During  Reverse  Transformation  (£  < 0). 

Next,  consider  the  case  of  reverse  transformation.  Using  (4.73), 

d$r 


da  = D : {de  — ocdT  + 


da 


(4.109) 


(4.110) 


(4.111) 
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The  consistency  condition  (4>  = 0)  is  obtained  by  differentiation  of  the  transformation 
function  such  that 


<9$  9$  <9$ 

dty*,T,s)  = — -.d<T+—<tr+—dt  = o 

da  dT  <9£ 

By  substituting  d a from  (4.111)  into  (4. 1 1 2),  i.e., 

] <9$ 

— :D:he-adT+w-  dA  + wdT+-dt  = 0 


(4.112) 


(4.113) 


: D : de  + < — — — — — : D : a > dT 


da 


dT  da 


, d$  ^ d$  d$)  „ n 
+ ^ — — : D : — 1 — — — / d£  — 0 


da 


da  <9£ 


(4.114) 


Arranging  for  d£  results  in 

<9<f> 


d£  = 


da 


\ <9$  d<f>  ) 

:D:rff+W_^:D:arT 


(4.115) 


where 


L f ^ <9$  <9$ 

b :=  < tt—  : D : — + 


(4.116) 


v<9cr  da  <9£  J 
Now  (4.1 15)  can  be  used  to  eliminate  <i£  in  (4.1 1 1)  and  obtain  the  relationship  between 
stress,  strain  and  temperature  increments  as  follows: 


da  = D : l de  — a dT  + 


<9<f> 

da 


d£. 


= D : d e — D : a:  dT  + D : 


da 


dT, 


= D : d e — D : a dT 


-!(»  £)(£  » '-)U 


f / (d®  d®  _ \)1„ 

“ I (D:  d^J  [df~  ~da  : :Q)|  bdT 
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= D : de 


: D : de 


1 

b 


D : cxdT 

( \ 1 


da  dTb 


dT 


+ D: 


<9$ 


da  / \ da 


= D 


D : 


- D-  D: 


da  ) 
d<$> 


D : aj 
( D :9i 


da 


D 


da  ) 
3$ 


-dT 

)sl 


d e 


da 


a dT 


- D 


— - -dT 


= D-  D 


da  J dTb 
<9<J> 


D - D : 


da  J 
d 4> 


D 


da 

C : de  + 0 dT 


D : 


<9$  \ 

da  J 

d<P 


- > : d e 


1 


da  I b 


: a + D : 


a «}- 


(4.117) 


The  continuum  tangent  material  modulus  tensor,  Cr,  and  continuum  tangent  thermal 
modulus  tensor,  0 r,  for  reverse  transformation  are  then  given  by 


C r = D - B <8)  B 


(4.118) 


3$  1 

@r  = ~£r ' a ~ &rt/bB 


(4.119) 


where 


, f 3$  _ 3$  3$ 

:d~a+~d^ 


„ / 3$\  1 

and 


(4.120) 
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4.3  Discrete  Time  Constitutive  Model  and  Implicit  Backward  Euler  Rule 

The  numerical  solution  of  nonlinear  boundary-value  problems  in  solids  mechan- 
ics is  based  on  an  iterative  solution  of  a discretized  time  version  of  the  corresponding 
boundary-value  problems.  Here  we  focus  on  only  the  constitutive  equations  to  be  com- 
puted and  implemented  into  the  element  stiffness  matrix.  The  numerical  algorithm  for 
integrating  the  rate  constitutive  equations  is  called  a stress  update  algorithm.  For  simplic- 
ity we  start  with  small  strain  hypothesis. 

4.3. 1 Basic  Algorithm  Setup  for  Strain-Driven  Problem 

Let  t be  a fictitious  time  quantity  such  that  t 6 [0,  T],  and  em  is  called  the  inelastic 
strain,  e.g.,  plastic  strain  or  phase  transformation  strain  in  SMAs.  The  real  relationship 
between  and  e™  is  as  follows: 


This  method  has  an  advantage  because  en+i  depends  only  on  known  values.  But  conver- 
gence rate  in  numerical  iteration  schemes  is  usually  slower  than  implicit  method.  The 
implicit  method  is  conceptually  more  consistent  because  the  current  strain  depends  on  the 
current-strain  rate  (and  previous  strain): 


(4.121) 


For  a numerical  solution,  an  approximate  relationship  between  em  and  em  is  required. 
The  usual  approximation  is 


(4.122) 


The  problem  is  that  em  is  changing  in  a time  interval,  and  we  do  not  know  which  value 
of  em  to  use:  at  the  start  or  end  of  the  time  interval,  or  some  intermediate  value.  The 
simplest  discretization  relies  only  on  known  quantities,  and  is  thus  explicit 


(4.123) 


(4.124) 
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This  equation  at  each  time  will  be  more  difficult  to  solve  because  and  e^n+1  are  both 
unknowns.  A more  general  method  can  be  constructed  that  incorporates  aspects  of  both 
the  explicit  and  implicit  methods  as  follows: 

e»+i  = en”  + [1  _ UJ\e™/\t  + ue™+1At  (4.125) 

where  the  explicit  method  can  be  obtained  by  setting  u — 0 and  the  implicit  rule  by 
u>  — We  call  the  semi-implicit  method  when  uj  — 1/2. 

For  the  system  of  ordinary  differential  equations  given  by  (4.72)  for  the  evolution 
of  phase  transformation  strain 

9$  • 

tP=JT~  £ (4.126) 

OCT 

we  can  obtain  the  evolution  equations  for  the  transformation  strain  which  are  discretized 
using  the  generalized  trapezoidal  rule  as  (Simo  and  Hughes  [1998,  section  3.6]): 

*?«=  6?  + [f„+i-f„]{[l-u.]^+a,^i]  (4.127) 

where  uj  = [0, 1]  and  subscript  n indicates  function  evaluation  at  time  tn  (similarly  £„+1), 
assuming  time  increment  At  = £n+1  — tn.  For  simple  choice  of  uj,  we  have  two  types  of 
integration  rules: 

• uj  = 0,  forward  (explicit)  Euler  integration  rule  (FER): 

epn+ 1=  e?  + [£n+i-£n]f^  (4.128) 

O <T 

• u — 1,  backward  (implicit)  Euler  integration  rule  (BER): 

e £+1  = 6 Pn  + [£n+l  - £n]  (4. 1 29) 

O (T 

A return-mapping  algorithm  consists  of  two  steps:  in  the  first  part,  a purely  elastic 
trial  state  (elastic-predictor)  is  computed,  usually  involving  an  excursion  away  from  the 
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yield  surface;  in  the  second  part  if  the  trial  state  violates  the  consistency  condition,  then 
an  inelastic-correction  is  computed  (inelastic-corrector  or  transformation-corrector)  us- 
ing the  trial  state  state  as  an  initial  condition.  During  the  process  of  the  time-continuous 
model  we  assumed  the  stress  as  control  variable,  whereas  during  the  process  of  the  time- 
discrete  model  we  assume  the  strain  as  control  variable  so  that  the  time-discrete  problem 
is  considered  to  be  a strain-driven  problem  for  which  a return-mapping  algorithm  is  used 
as  the  solution  algorithm  (Simo  and  Hughes  [1998]).  In  other  words,  given  the  incre- 
ments of  strain  field,  {Ae„+1}  at  time  tn+1,  the  task  is  to  update  the  field  variables  from 


In  analogy  to  the  case  of  plasticity,  the  trial  state  is  in  general  defined  as  the  state 
obtained  from  the  previous  solution  state  (cr„,  e„)  taking  into  account  the  new  strain 
( e n+i)  and  assuming  no  evolution  for  the  phase  transformation  quantities.  That  is,  when- 


transformation. The  enforcement  of  the  consistency  condition  makes  it  possible  to  com- 
pute the  admissible  solution  and  the  evolution  of  the  phase  transformation  quantities. 


ever  the  trial  state  ( cr  ^x)  violates  the  consistency  condition,  the  step  is  on  the  phase 


(4.130) 


(4.131) 


where 


(4.132) 


The  above  equations  can  also  be  rewritten  as  follows: 
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- atry 

— ” n+ 1 


= 


ftn+l 

/ D :dept 
Jtn 

3$ 


da 


n+1 


da 


(4.133) 


where  <&(  a , T,  £)  is  referred  to  the  transformation  function.  Geometrical  illustration  of 
this  return-mapping  algorithm  is  presented  in  Figure  4.5. 


4.3.2  Implicit  Return  Mapping  Algorithm  of  Stress  Integration 

In  order  to  integrate  the  transformation  flow  rule  for  the  transformation  correction, 
we  use  backward  Euler  method  ( uj  = 1)  resulting  in  a nonlinear  algebraic  set  of  equations, 
which  are  solved  using  Newton’s  iteration  method  (Simo  and  Hughes  [1998,  p.  143]).  We 
thus  have  a set  of  constitutive  equations  as  follows: 

£?+!=£?  + [fn+1  - £n] (4. 1 34) 

o a 

and 


<H,  = D : {<r„+1  - a [T„+i  - To]  - «*£»}  (4.135) 

constrained  by  the  discrete  Kuhn-Tucker  conditions. 

For  the  forward  transformation  (A  — > M),  the  discrete  Kuhn-Tucker  conditions  are 
expressed 


[£n+l  £n]  > 0, 


& n+l)  d^n+l-i  £n+l)  5:  0) 

[fn+i  - &]  •$/(<rn+1,TB+1,U)  = 0 (4.136) 

and  the  discrete  consistency  condition  using  A£n+i  :=  £n+1  — £n  = 0 is 

^/(^n+1,Tn+1,en+l)=0  (4.137) 

For  reverse  transformation  (M  —>  A),  the  discrete  Kuhn-Tucker  conditions  are 
expressed 


[Cn+l  £n]  ^ 0, 
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^*r(  & n+l)  ^ri+l?  £n+l)  A 0, 

[£ri+l  — £n]  *^*r(  & n+ 1>  Tn+ li  Cn+l)  = 0 (4. 138) 

and  the  discrete  consistency  condition,  using  A£n+1  :=  £n+1  — = 0,  is 

(I*r(  <j  n+i,  £n+i ) = 0 (4.139) 

At  the  A>th  iteration  step,  if  the  transformation  criterion  is  satisfied,  i.e., 

if  :=  $(  <r  S1?  Tn+1,  e„)  < 0 (4. 140) 

then  the  trial  thermoelastic  state  is  the  final  state.  However, 

if  ^1:=$(<r?]1,Tn+1)e„)>0  (4.141) 

then  the  Kuhn-Tucker  conditions  are  violated  and  the  trial  state  lies  outside  the  trans- 
formation surface.  Transformation  correction  employing  backward  Euler  rule  (BER)  for 
flow  rule  is  then  used  to  restore  consistency  and  the  solution  j a ^°|l5  Tn+1,  £n  j are  taken 
as  the  initial  conditions  for  the  transformation  corrector  step. 

4.3.3  Elastic  Prediction  and  Transformation  Correction 

Physically  the  trial  state  is  obtained  by  freezing  the  transformation  flow  due  to  the 
transformation  during  the  (n  + l)-th  loading  incremental  step,  i.e.,  Afn+1  = 0.  The 
thermoelastic  prediction  problem  is  given  by 

e n+\  = e„  + Aen+1  = in  en+(i0)  = e?  (4.142) 

Thus,  the  trial  stress  at  k = 0 can  be  expressed 

<4+1  = D : {€„+,  - a[r„+,  -71]  - ej?}  (4.143) 

For  k > 0,  i.e.,  from  the  second  iteration  step, 

o-nh  = D : { £„+1  - a [T„+1  - T„]  - £*<,*>}  (4.144) 
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If  ^L+i  > 0,  the  transformation  correction  procedure  is  adopted  at  the  (n  + 1)- 
th  loading  increment  for  (4.134)  and  (4.135).  It  is  important  to  note  here  that  replacing 

(k) 

cr  n by  a implies  that  the  implicit  integration  of  the  transformation  strain  flow  rule  is 
carried  from  the  initial  thermoelastic  prediction  to  the  final  state. 

The  nonlinear  system  of  these  algebraic  equations  is  valid  by  defining  the  trans- 
formation condition  valid  for  the  transformation  corrector  step  and  transformation  strain 
residual,  R%+i\  based  on  implicit  backward  Euler  rule  (BER)  from  (4.134)  for  the  fc-th 
iteration  as  follows: 


:=  4>(^Si .r„+i,0 


(fc) 


:(*) 


(4.145) 


R 


Pt,(k) 

n+1 


rPt’ik)  , ,pt 

■—  en+l  ' en 


!,,(*>  _ , , 

q LSn+1  snj 


(4.146) 


We  want  and  R^+i1  to  be  zero  at  the  end  of  the  iteration  process  (Simo  and  Hughes 
[1998,  p.144];  Belytschko,  Liu  and  Moran  [2000,  p.280]).  That  is, 


,pt,(k) 


lim  $ 

k— »oo 


(k) 

n+1 


0 and 


lim  0 

k— >oo 


(4.147) 


Linearization  of  (4.145)  and  (4.146)  now  leads  to  the  following  way 


+ 


da  'Act”+>  + 


(*) 


dt 


a eSi  = o 


(4.148) 


T}Pt,(k)  A cPt’ik)  , 

n+1  U£«+l  “r 


d cr 


(*) 

n+1  a c(k) 
^?n+l 


+ 


[d+1  - Zn}^^-  ■ =0,  (Z>  0) 

o cr  L 


(4.149) 


0,  (Z  < 0) 

where  (4.149)  is  obtained  by  the  Kuhn-Tucker  condition  when  — £n]  ^ 0.  The 


two  equations  (4.148)  and  (4.149)  provide  7 equations  to  solve  for  13  unknowns,  i.e., 

£6x1  and  Z € Rlxl. 


{Ao-Jgj,  Ae^1,  A^}.  Note  that  <r  G M6xl,  e*  € 
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Remark  4.6.  Note  that  = 0 for  reverse  transformation  (M  — > >1)  because 


a<t> 


<*) 

n+l 


during  the  reverse  transformation  remains  constant  with  its  components  which  are 
dependent  on  the  transformation  strain  at  the  end  of  the  previous  forward  transformation. 
Recall  that  for  the  case  of  J2  flow 

f| zl\ct\  (£  > 0) 


dcr 


-pt  £ 


i ept 


■pt  c r , 


(e<o) 


(4.150) 


The  rest  of  the  6 equations  are  obtained  by  taking  the  increment  of  the  stress  strain 
relation  (4.135): 


o’n+i  = D 1 : 


a[Tn+1-T0]~  ef 


;"+l  ^ L^n+l  “ ^OJ  - cn+1 

Expressing  it  for  the  total  increment  of  stress,  we  obtain 


(4.151) 


Acrn+1  = D 


A e n+i  — ot  A Tn+i  — A 


: Pt 
■ n+1 


(4.152) 


Using  above  equation  and  noting  that  e and  T are  fixed,  i.e.,  A e = AT  = 0,  during 
the  return-mapping  (transformation  corrector)  step,  (4.152)  can  be  written  in  incremental 
form  for  A e (f+1  at  the  k- th  iteration  step  as  follows: 


Ae^,  = — D_1  : A cr[k) 


' /,n+ 1 


n+l 


(4.153) 


Three  equations  (4.148),  (4.149),  and  (4.153)  form  of  a system  of  equations  to  solve  for 

{A  cr  Ae^f,  A^jj  for  both  forward  and  reverse  transformations. 

During  Forward  Transformation  (£  > 0). 

During  forward  transformation  (A  ->  M),  re-arranging  (4.149)  for  Aejf  gives 

>(*) 

■ n+ 

dV 

and  subtracting  from  (4.153)  yields 


\ pt,(k)  _ j3Pt,(k)  dtf?n+1  Af(k)  rc(fc)  c A 

LXef,n+ 1 - "n+l  + A4n+ 1 + (£n+l  ~ £nj  g 2 • A, 


.(*) 

n+l 


(4.154) 


^ " /,n+ 1 — “n+l 


n+l 


3a  A^n+1 


(4.155) 
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in  case  of  J2  flow 


1 

6»cr 


= +A 


(k) 
n+ 1 


(4.156) 


and  the  algorithmic  elastic  stiffness  tensor  is  defined  as 


cl*)  .. 

^n+1  •" 


d-1  + [d‘A  - 

OCT  Z 


-1  -1 


(4.157) 


This  expression  is  called  exact  Hessian  matrix  (Simo  and  Hughes  [1998,  p.  144]).  Gradi- 
ent of  transformation  direction  tensor  with  respect  to  the  stress  is  given14 


<924> 

5^ 


12  6l_ 

3 llsll 


I 


-^1  ® 1 — Tj— 7 

3 s 


s 

Isl 


(4.158) 


where  eL  is  maximum  uniaxial  transformation  strain,  I is  the  4th-order  symmetric  identity 
tensor,  and  1 is  the  second-order  identity  tensor. 

Substituting  (4.155)  into  (4.148),  the  increment  of  the  martensite  volume  fraction 
is  obtained  for  the  A;-th  iteration  as 


$ 


(*) 
n+ 1 


+ 


<9$ 


(k) 

n+ 1 


da 


w(*) 
,Jn+ 1 


R 


pt,{k) 

n+1 


~~da 


A c(k) 
^Sn+l 


(k) 

+ = 0 


' Tfpt,(k ) _ — .(fc) 

^n+1  • -rtn+l  1—1 n+1 


da 


(k) 

"+1Ad+i 


°^n+l  A/-(fc)  _ n 
— u 


4> 


(fc) 

n+1 


d$(k) 


n+1  _ .^(fc) 


da 


,pt,(k) 


. R 1 

n+1  • ■rLn+l 


^n+l  -.(fc) 

da  ' “”+1 


5$n+l  Ac(k)  , d^nh  Af(k)  n 

da  ^+i  + -^-^n+1-0 


4> 


(fc) 

n+1 


d <J> 


(fc) 

n+1  . — .(&) 


da 


. nPt,(k) 
n+1  • ■r*'n+l 


Sfgi  =w 

da  ‘ “n+1 


<94> 


(fc) 

n+1 


3*iSi 


ao- 


Ad+1  = 0 


A ,c(fe)  _ J adfc)  ^/.n+1  . — ,(fc)  . jypt,{k)  l 1 

‘-iS/,n+l  _ 4 Vfn4.1  • — n-Ll  • > 77" 


■ /,n+l 


da 


■‘n+l 


' n+1 


(fc) 
ln+ 1 


(4.159) 


14 


Derivation  will  be  provided  in  Appendix. 
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where13 


5(*)  ._ 
an+ 1 •— 


(k) 

/,n+ 1 . — i(fc) 

da  ' “"+1  ' 


Oa 


(4.160) 


During  Reverse  Transformation  (|  < 0). 

A similar  analysis  can  be  performed  during  reverse  transformation  (M 
arranging  (4. 149)  for  A e 


; n+l 


A cf*1’!1)  _ T>Pt<W  I 
^ r,n+l  -*^77.4-1  i 


n+l 


(k) 

n+l  A 


<9<t 


Ae 


n+l 


A).  Re- 


(4.161) 


and  subtracting  (4.161)  from  (4.153)  yields 


Ao-«  i = D 


dk) 


r>  Pt,(k)  , ^V.n+l  A c(k) 
ltn+ ! ' ^U+l 


(4.162) 


in  case  of  J2  flow 


apr,n+ 1 

da 


( k ) 

■ n+l 


(4.163) 


Substituting  (4.162)  into  (4.148),  the  increment  of  the  martensite  volume  fraction 
is  obtained  for  the  k- th  iteration  as 


+ d®n+l 
n+l  ' 

da 


D 


- R 


Pt,(k)  d$l+i  AAk) 

n+l  ' a ^?n+l 

cl  a 


= 0 (4.164) 


Cl  + 


Cl 

da 


D : R 


Pt,(k) 

n+l 


da 


D : ^±i A^i  | + ^|±I  A^x  = 0 


*«  - 
n+l 


(k) 

71+1  ■ D : Rpn^k)  + 


d$ 


(k) 

n+l 


3®1S,  ,M) 


3$ 


(*) 


da 


di 


: D : AAEtiACi  + A^aCi  = o 


<9cr 


rf>(fc)  °^n+l  . p,  . TfPt,(k) 

^n+1  U Un+1 


d+k) 


n+l 


da 


d<b^ 

pv  . u'*n+ 1 

' da 


+ 


3$ 


(k) 

n+l 


ae 


Ad+i  = 0 


15  See  Simo  and  Hughes  [1998,  p.  145]  for  plasticity  problem.  In  fact,  the  above  expression  is  the  same 
as  the  plasticity  problem  case. 
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Solving  for  A 


where 


(*)  _ 
/,n+l 


$ 


(*) 

r,n+l 


^^r.n+l 

<9(7 


: D : i? 


(4.165) 


/9(f)  (fc)  /J<T,(fc)  /9<T>(fc) 

^T-.n+l  . j-v  . °^r,n+ 1 , g^V,7i+l 

da  da  <9£ 


(4.166) 


Update  procedure. 

In  (4.153),  (4.155),  and  (4.159)  for  forward  transformation  and  (4.153),  (4.162), 
and  (4.165)  for  reverse  transformation,  we  have  a complete  set  of  equations  for  the  13 
unknown  variables  {Aer^,  Ae£+*\  Af^+jj.  The  next  step  is  to  update  the  stress, 
transformation  strain,  and  martensite  volume  fraction  by 


,(fc+i) 

71+1 


.pt,(fc  + l) 
' 71+1 


£.(fc+l)  _ 
Sti+1 


(*)  + A rr  {k) 

(4.167) 

pt,(fc)  A pt,(k) 

e7i+l  _r^ie7i+l 

(4.168) 

d+i  + A 

(4.169) 

4.3.4  Geometric  Interpretation  of  Closest-point  Projection  Algorithm 

Here  we  examine  the  geometric  interpretation  of  the  closest-point  projection  method 
(CPPM)  associated  with  the  discrete  constitutive  relations  (Simo  and  Hughes  [ 1 998,  p.39, 
118]).  Substituting  (4.134)  into  (4.135),  we  obtain  for  the  (k  + l)-th  iteration  at  tn+x 

,d$ik+1) 


<’•&*’  = D : { e„+1  - a [T„+1  - To]  - 


Recall  (4.143): 


(4.170) 


.(0) 

71+1 


D : { en+1  — a[Tn+1-T0]-  e*} 


(4.171) 


Subtracting  the  resulting  equation  yields 


D-:(, 


. (k+1) 

71+ 1 


fy  (°)  ) _ _rc(*+D  _ c 

17  71  + 1 ) — [Sti+I  S7i 


(fc+1) 

71+1 


da 


(4.172) 
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Figure  4.5.  Geometrical  illustration  of  closest-point  projection  . 


Finally, 


(4.173) 


The  return  (or  relaxation)  path  of  the  stress  from  the  elastic  prediction  with  each  iteration 
is  not  arbitrary  but  is  defined  by  the  projection  of  the  thermoelastic  prediction  on  the 


procedure  in  such  a manner  that  at  (k  + l)-th  iteration,  in  the  six-dimensional  stress  space 


4.3.5  Algorithmic  (consistent)  Tangent  Moduli 

The  important  next  step  is  that  the  tangent  moduli  appearing  in  the  linearized  prob- 
lem due  to  Newton-Raphson  method  must  be  obtained.  In  continuous  time  problem,  the 
continuum  tangent  material  modulus  (CMAM)  and  continuum  tangent  thermal  modulus 
(CTEM)  result  from  the  enforcement  of  the  classical  consistency  condition  on  the  con- 
tinuum problem.  In  order  to  maintain  the  quadratic  rate  of  asymptotic  convergence  in 
Newton-Raphson  method,  the  notion  of  algorithmic  (or  consistent)  tangent  material  mod- 
ulus (AMAM)  and  algorithmic  (or  consistent)  tangent  thermal  modulus  (ATEM)  arises 


transformation  surface,  That  is,  the  stress  state  is  updated  during  the  iterative 


defined  by  the  metric  of  elasticities  D,  the  stress  state  er  lies  on  and  the  vector 

joining  er  ■ That  is,  the  stress  will  be  updated  normal  to  each  iterative  transformation 
surface  in  stress-space.  In  other  words,  the  vector  joining  er  ^ to  er^  is  normal  to 
the  transformation  surface  (Figure  4.5). 
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due  to  the  enforcement  of  the  consistency  condition  on  the  discrete  algorithmic  problem 
(Simo  and  Taylor  [1985];  Zienkiewicz  and  Taylor  [1991]). 

Recalling  the  constitutive  relation  from  the  time  continuous  model, 

cr  :=  p—  = D : { e - a [T  - T„]  - ept } (4.174) 

The  discrete  model  of  stress-strain  relations  are  from  (4.134)  and  (4.135)) 

ern+1  = D : {e„+1  — a[Tn+1  -T0]  - ept+1} 

= D : { e n+1  - a [Tn+1  - T0]}  - D : e jf+1  (4.175) 

where 


en+i  = ep„‘  + [fn+i  - f„]  An+1  (4.176) 

Here  can  be  obtained  by  using  the  discrete  transformation  function  $n+1. 
Algorithmic  Tangent  Material  Modulus  in  FPT. 

First,  we  linearize  the  stress-strain  equations  by  rewriting  it  at  time  n + 1 by  taking 
the  total  differential  of  the  discrete  flow  rule  and  the  constitutive  stress-strain  relation  as 
follows: 

We  now  linearize  the  equation  of  stress-strain 


dcrn+l  = D 


d £ n+i  <*dTn+1  ds1 j_|_j 


(4.177) 


Using  the  linearized  form  of  (4.134)  for  forward  transformation 


pt  d$n+1  . r„  _ , <924> 

d € „ , i — 


• n+1  a 

da 

Arranging  for  d cr  n+1  gives 

»-i  , ^ r id2$n+ 1 

da 2 


d^n+l  T [£n+l 


n+1 


da 2 


da  n+i  (f>0)  (4.178) 


D + [£n+l  — £n 


d$f  4.1 

: dan+l  = den+1-  cxdTn+l J^±±d^n+1(4.n9) 

da 
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For  forward  transformation  (A  —>  M),  an  algorithmic  relation  for  do  n+l  in  terms  of 
d e n+1,  dTn+i,  and  d£,n+ x is  obtained  as  follows: 


d o n+i  — “ 


n+l  • 


d e 


n+l 


a dT, 


n+l 


/,n+ 1 


da 


df 


n+l 


(4.180) 


where  the  algorithmic  elastic  stiffness  tensor  is  defined  as 

,d2$n+l 
do 2 


“n+1  • — 


D + [£n+l  — £n]‘ 


(4.181) 


Next,  in  order  to  determine  dfn+i  > 0,  we  want  to  determine  the  discrete  trans- 
formation multiplier,  <7£n+i  ■ Differentiation  of  the  discrete  consistency  condition  (4.137) 
yields 


0<l>t 


d$f,n+ 1 ( <r  , T,  0 = " : do  n+1  + ^V+i  dTn+1  + d£n+1  = 0(4. 1 82) 


do  qt 

Now  substitution  of  do  n+i  in  (4.180)  into  (4.182)  yields 


d®n+ 1 , d$n+1  , d$n+ljr 

■ d o n+i  -| ^—dln+i  -| — — dfn+ 1 = 0 


do 


dT 


df 


(4.183) 


n+l 


do 

5$n+ 1 


Jn+ 1 


, <9$n+l 

d c n+i  otdln+ 1 Q _ aqn-)_x 


dTr 


5$. 


n+l 


n+1 


^Cn+1  — 0 


(4.184) 


<9$n+1  r 

^ | “n+1  • u 6 n+1  “n+l  • ^ dTr i _i_ i — 

5$n+1 


<9$ 


n+l  “n+1 


n+1 


dcr 


n+l 


H -ffff-dTn+i  + ^%±l^B+i  = o 


dT 


dl; 


(4.185) 


3$n+l  . „ , 5$n+ 1 „ 

“n+i  • d e n-j_i 


do 

d$n+ 1 

dT 


t/T, 


n+l 


n+1 


d o 


dfn+l  — 0 


n+1 


. ci  dTn+i 


<9$ 


n+1  , 


<9$ 


n+l 


do 


“n+l 


do 


df 


n+l 


(4.186) 
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3$n+i  ^ j d$n+l  jrr  <9$ 

• “77+I  • aen+ 1 i ctm  aln+l 


77+I 


da 

'd$n+1 


da  ‘“n+1'  da 

and  arranging  for  d£„+1  results  in 


dT 

<93>n+i  d$n+1 


d£ 


da 

d^n+l  ~ 0 


“77+1  ■ 


ct  dT, 


71+1 


(4.187) 


d£n+l  — 


^77+1.  „ , 

Q & ' 1— '77+1  • d & 77+1  I 


■ “77+1  • d 6 n+l 


<94> 


77+1 


dT 


d<f>n+l  . „ 

a • “77+1  ■ C* 

da 


dT 


77+1 


*77+1 


+ 


da 

^77+1 

dT 


®77  + l 


d$n+ 1 . „ 

• “77+1  • 

da 


®77+l 


-dT 


77+1 


(4.188) 


where 


®77  + l • 


d<S> 


/,77+l  _ , 


<9$ 


da  ‘~n+1’  c)ff 
Consequently  substitution  of  (4.188)  into  (4.180)  yields  the  expression 

<937,77+1 


/,77+l 


<9$ 


/,77+l 


(4.189) 


d & “ *77  + 1 • d 6 77^-1  “77  + I . a dT 


77+1  “77+1 


<9cr 


d£, 


77+1 


(4.190) 


“77+1  ' ds  7,4-1  £!■ 


77+1  “77+1  • 


ct  dT7 


77+1 


“77+1  • 


77+1 


7>b 


77+1 


+ 3 


77+1 


da 

^77+1  ' 


da 


■ “77+1  ■ d 6 


1 


77+1 


'3$ 


da  ) 

“77+1  • de  7J  + 1 — 

r 7 d$n+1 


77+1  _ 

“77+1  : a 


®77  + l 

5$„+l  \ 1 


da 

d^+i 

da 


dT  ) an+i 

d^n+i  „ 


dT, 


<9cr 


■*77+1 


: de 


+ U“”+1'  Off 

“77+1  ■ dT) 7+1 

1 “77+1  ( “77+1 


<93>77  + l „ 

— : £177+1  : ol 

da 


d<l> 


77+1 


77+1 


1 


77+1 


®77  + l 


dT 


77+1 


da 


“77+1  • 


~di 


77+1 


®77+l 


®77  + l 


: de 


dTr 


77+1 
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“n+l 


0$n+l\ 

da  J 


d$n+l\  1 | 

da  J an+i  j 


■ ^ dTn_i_i 


(~  ^n+l^  <9$n+l  1 

rn+1  • <9<r  J dT  an+i 


dTn+\ 


(4.191) 


Defining 


£ 


da 


f,n+ 1 • — 


de 


and  © 


da 


n+l 


f,n+ 1 ■ — 


dT 


n+l 


(4.192) 


the  algorithmic  tangent  material  modulus  (AMAM),  Cf,  and  the  algorithmic  tangent  ther- 
mal modulus  (ATEM),  0 f,  for  the  forward  phase  transformation  ( A — ► M)  are  derived 
as 


£/,n+ 1 “n+l  A 


n+l 


Ln+1 


(4.193) 


where 


© /,n+l  — Tj 


<9$, 


®n+ 1 • 


'<9$ 


/,n+l 


da 


A-n+1  • — ( “n+l 


dT 

/= An+1 

V®n+ 1 

<9$/,  n+l  <9<f>/„+ 1 

"1  * a 

OCT 

U s 

' CQ 

1 

1 

d&f,n+ 1 \ 

1 

da  ) 

y/O-n+1 

(4.194) 


(4.195) 


(4.196) 


k— 'n+l  • 


D'1  + [f, 


,g2i 


n+l 


n+l 


da 


(4.197) 


It  is  observed  that  as  the  time  step  A tn+1  :=  tn+1  - tn  -*•  0 and  A£„+1  :=  £n+1  - -> 

0,  i.e.,  the  consistent(or  algorithmic)  tangent  stiffness  tensors  reduce  to  their  continuum 
counterparts  as  At„+1  -*•  0 in  (4.1 10)  (Simo  and  Hughes  [1998];  Qidwai  and  Lagoudas 
[2000a]).  The  symbol  tilde  over  the  letter  means  the  quantity  in  algorithmic  part. 
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Algorithmic  Tangent  Material  Modulus  in  RPT. 

First,  we  linearize  the  stress-strain  equations  by  rewriting  it  at  time  n + 1 by  taking 
the  total  differential  of  the  discrete  flow  rule  and  the  constitutive  stress-strain  relation. 


d a 


n+1 


= D 


d e n+i  — at  dTr 


n+1 


d e 


pt 

71+1 


(4.198) 


Using  the  linearized  form  of  (4. 134)  for  reverse  transformation 

de^+1  = -^d^n+lt  (i<  0) 

OCT 


(4.199) 


Arranging  for  d a n+1  gives,  for  reverse  transformation  (M  ->  A),  an  algorithmic  relation 
for  d a n+1  in  terms  of  d e n+1,  dTn+i,  and  d£n+i  is  obtained  as  follows: 


d cr  n+i  — D : j d e „+i  — a dTn+\  H ^£n+i 


(4.200) 


Next,  we  want  to  determine  the  discrete  transformation  multiplier,  d£„+ 1.  Differen- 
tiation of  the  discrete  consistency  condition  (4. 1 39)  yields 


^7/f,  / T c\  ^^r.n+l  ; _ , d^r.n+l 

d^r,n+ 1(  & > a j £)  — — q : d CT  n+1  + 


do-  dr 

Now  substitution  of  da  n+1  in  (4.200)  into  (4.201)  yields 

d$n+i  + _ , d®n+ijrr  , d$n+ \Jt. 

• da  n+i  H ^ dl  n+i  H — — a£n+i  = 0 
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d$n+i  . ~ , <9$n+i 
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and  arranging  for  d£n+1  results  in 
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Consequently  substitution  of  (4.207)  into  (4.200)  yields  the  expression 
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the  algorithmic  tangent  material  moduli  tensor  (AMAM),  £r,  and  the  algorithmic  tangent 
thermal  moduli  tensor  (ATEM),  © r,  for  reverse  transformation  ( M — > A),  are  derived  as 
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’n+l 


It  is  observed  that  for  the  reverse  transformation,  the  consistent  (or  algorithmic)  tangent 
stiffness  tensors  are  the  same  as  their  continuum  counterparts  in  (4.120). 


4.4  Determination  of  Material  Parameters 


Descriptions  on  input  data  in  this  section  are  based  on  Funakubo  [ 1 987]  and  Lagoudas 
and  Bhattacharya  [1998],  The  proposed  transformation  hardening  function  is  of  the  su- 
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perposition  of  linear  and  Voce-type  exponential  hardening  (recall  from  section  4.2.6). 


1 


/f(£)  ==  --(pAr,0)(Ms-Mf)e 


- (pAVo)  (Ms  - Mf ) 

+mi( 


e+Vexpt-ZJ"?) 
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//?(  0 := 
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(pA77o)(A/  - As)£2 


(pArjo)  (Af  - As) 


-m  if 


f + ^expH^f) 
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and  the  derivative  of  hardening  function  with  respect  to  f is 

0/fc"(O 


(4.215) 


(4.216) 


= —(pAp0)(Ms  — Mf)£ 

- (pAr]0)  (Ms  - Mf)  (l  - exp(-/?A/f))  + (4.217) 

= —(pAr]o)(As  - Af)Z 

~ {pAp o)  (A/  - As)  (l  - exp(-^O)  - mi  (4.218) 


which  provides  both  linear  and  nonlinear  isotropic  hardening  term  characterized  by  the 
isotropic  hardening  exponent  / 3M  and  (3A. 

The  maximum  uniaxial  phase  transformation  strain  eL  can  be  obtained  from  a uni- 
axial test  by  measuring  the  residual  strain  retained  after  the  specimen  has  been  fully 
loaded  to  martensite  phase  followed  by  complete  unloading.  For  1-D  case,  A = eL.  The 
slope  of  the  curves  in  stress-temperature  phase  diagram  during  forward  transformation 
and  reverse  transformation,  respectively,  can  be  observed  phenomenologically. 

- -(sr 
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Using  the  above  given  material  parameters,  we  can  determine  the  rest  of  material 
material  parameters  in  the  following  steps: 

• In  order  to  determine  the  difference  in  entropy  pArj0,  recall  from  the  previous  sec- 
tion the  transformation  function,  similar  to  yield  function  in  plasticity 


:=  aeL  + (pAr)0)T  - (pAE0)  - (-pAr]o)(Ma- Mf^-rri! 
+ (pAVo)  (M,  - Mf)  (l  - exp(-/?M0)  - Yf*  (4.220) 

:=  aeL  + (pAVo)T-(pA£o)-(-pAVo)(Af-As^  + m1 
+ o)  {Af  - As)  (l  - exp(-/^f))  + Y*  (4.221) 


For  the  case  of  1-D  uniaxial  tension,  the  stress-temperature  phase  diagram  based  on 


ure  4.6  (Brinson  and  Lammering  [1993],  Ford  and  White  [1996],  Lagoudas  et  al. 
[1996],  Auricchio  [1995]).  The  entropy  difference  pAr]0  per  unit  volume  between 
the  phases  can  then  be  determined  phenomenologically  by  the  slope  (^)  of  these 
curves  at  zero  stress.  To  analytically  obtain  the  slopes,  we  first  take  the  differen- 
tial of  (4.220)  for  constant  martensitic  volume  fraction  £ along  the  transformation 
surfaces,  the  slope  of  the  curves  can  be  obtained  by  the  following  equation: 


(?  >0) 

(I  < 0) 


(4.219) 


We  consider  1-d  uniaxial  tensile  case.  Then 


experimental  data  for  fixed  martensitic  volume  fraction  (£  = 0)  are  plotted  in  Fig- 


0 


(4.222) 
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Figure  4.6.  Stress-temperature  (cr  — T ) phase  diagram  of  superelasticity. 


then  it  becomes 

da_  = _J|  = (pAr]0) 

dT  f 

Then 


(4.223) 


(pA?7o) 


(4.224) 


The  above  equation  is  usually  called  the  Clausius-Clapeyron  equation  for  the  first- 
order  phase  transformations  under  uniaxial  stress.  The  variation  in  stress  needed  to 
produce  the  stress-induced  martensite  phase  transformation  increases  linearly  with 
temperature. 


Remark  4.7.  See  also  Funakubo  [ 1 987,  p.46]  and  Otsuka  and  Shimizu  [ 1 986,  eq.  1 9] 

for  detailed  description  about  Clausius-Clapeyron  equation.  If  the  pressure  in  the 

conventional  form  used  in  thermodynamics  is  replaced  by  the  uniaxial  stress,  then 

the  right  side  is  changed  from  positive  to  negative.  For  general  case  of  first-order 

phase  transformation,  the  Clausius-Clapeyron ’s  equation  is  defined 

dP  dq  A7712 

dT  dv  TAV 


(4.225) 
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where  A7 712  is  latent  heat  of  transformation  from  phase  1 to  2.  The  Clausius- 
Clapeyron  equation  relates  the  slope  along  the  coexistence  line  with  the  change  in 
entropy  and  volume  of  the  substance  as  it  crosses  the  line,  i.e.,  changes  phase.  I 

• The  parameters,  Y*,  ml5  and  pA£0  can  be  determined  after  considering  both  for- 
ward and  reverse  transformation  (Raniecki,  Lexcellent  and  Tanaka  [ 1 992],  Lagoudas 
and  Shu  [1999,  602]).  During  the  forward  transformation, 

<f>/(a,T,£)  :=  aeLY(pAp0)T  - (pA£0) -mx-Y* 

-(-pAp0)(Ms  - Mf)i  + (pAri0)  (Ms  - Mf)  (l  - exp(-/?M£)) 

= 0 (4.226) 

During  the  reverse  transformation, 

$r(a,T,€)  :=  aeL  + (pAp0)T  - (pA£0) + m1+Y* 

-(-pAp0)(Af  - Aa)£  + (pAp0)  ( Af  - As)  (l  - exp(-/?A£)) 

= 0 (4.227) 

Consider  the  specimen  with  setting  a = 0 is  being  cooled  down  from  full  austenitic 
phase  at  a high  temperature  to  the  temperature  Ms  and  then  Mf,  while  the  corre- 
sponding volume  fraction  is  zero.  First,  consider  the  case  of  forward  transformation 
(A  — > M ).  When  T = Ms,  £ = 0.  Thus 

{pApo)  Ms  - (pA£0)  - mi  - (1  - 1)  - F*  = 0 (4.228) 

and  when  T = Mf,  £ = 1, 

(pAp0)  Mf  - (pA£0)  - ( -pAp0)(Ms  - Mf)  - mi 

+ (pAp0)  (Ma  - Mf)  (l  - exp(-/?M))  -Y*  = 0 (4.229) 



nonlinear  term 
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Likewise  during  the  reverse  transformation  (M  — > A),  at  low  temperature  to  the 
temperature  As  and  Af,  the  corresponding  volume  fraction  of  martensite  £ reduces 
from  100%  to  0%.  When  T = As,  £ = 1.  Thus 


(pAp0)  Aa  - (pA£0)  - (—pArjo)(Af  - As)  + mx 
+ (pAp0)  (Af  - As)  (l  - exp(-/^))  +y*  = 0 (4.230) 

V v 

-v 

nonlinear  term 

and  when  T = Af,  £ = 0, 

(pAp0)  Af  - (pA£0)  + m!  - (1  - 1)  + y*  = 0 (4.231) 


• The  parameter  m i can  be  found  using  continuity  condition,  f{f(  1)  = f£(  1).  Thus 


/f(l) 
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1 + ^i7ex  P(-/^M) 
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(4.232) 


AA(i) 


1 
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(pArio)(Af  - As) 


~ (pApo)  (>1/  - A) 


—m  i 


(4.233) 


Thus, 

toi  = ^(-pApo)  [(A/  - Aa)  - (Ma  - Mf)]  + (rM  - r'4)  (4.234) 

• The  difference  of  internal  energy  pA£0  = p[8{f  - 8q]  where  8™  and  8*  are  the 
internal  energy  at  martensite  and  at  austenite,  respectively,  can  be  determined  as 
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follows  (Funakubo  [1987,  p.  1 3],  Raniecki  et  al.  [1992],  Lagoudas  and  Bhattacharya 
[1998,  p.360])  From  a purely  thermal  transformation,  i.e.,  a — 0 and  Y = 0,  the 
change  in  internal  energy  pA£0  can  be  determined 


Y(p,T, 


0 


0 

aeL  + (pAp0)  T - (pA£0)  - 
0 + 0 + (pApo)  T - (pA£0)  (4.235) 


which  can  be  simplified  as 


pA£0  = (pArjo)  Teq  (4.236) 

where  the  equilibrium  temperature  is  defined  as  Teq  = | [Ms  + A f], 16  Thus,  pA£0 
can  be 

pA£0  = ^ (pArjo)  [Ms  + A f]  (4.237) 


• The  internal  dissipation  energy  can  be  obtained  from  (4.230): 

Y*  = — (pApo)  As  + -(pArj0)(Ms  + Af)  — (pArj0)(Af  — As) 

-(pApo){Af  - As)  (l  - exp(-/5j4))  - mi  (4.238) 


16  In  Funakubo  [1987,  p.  1 3]  or  Berman  and  White  [1996],  T0  is  called  the  average  PT  temperature 
or  equilibrium  PT  temperature.  The  location  of  this  equilibrium  temperature  varies  widely  as  reported  in 
the  literature  (Sun  and  Hwang  [1993a];  Duerig  et  al.  [1990,  p.26];  Auricchio  et  al.  [1997];  Brinson  and 
Lammering  [ 1 993]). 


CHAPTER  5 

NUMERICAL  EXAMPLES 


In  this  chapter,  we  provide  numerical  examples  which  demonstrate  the  qualitative 
and  quantitative  behavior  of  the  numerical  implementation.  In  order  to  show  the  numer- 
ical efficiency  of  our  model,  the  constitutive  model  has  been  implemented  into  finite- 
element  formulation  of  low-order  solid-shell  element.  A set  of  algebraic  equations  result- 
ing from  finite  element  discretization  of  the  boundary  value  problem  are  solved  globally 
via  a Newton-Raphson  scheme.  Representative  simulations  are  presented  for  a solid-shell 
element:  Uniaxial  tensile  loading,  simply-supported  three-point  bending  problem,  and  a 
problem  of  asymmetry  effects  in  uniaxial  tension  and  compression.  Also  the  performance 
of  the  present  formulation  involving  large  deformation  in  shell  structure  is  also  tested.  All 
the  presented  examples  are  subjected  to  stress-induced  martensite  phase  transformation. 
Table  5.1  describes  material  parameter  inputs  which  are  used  in  computation.  Table  5.2 
and  5.3  show  material  properties  used  in  other  models  and  our  model.17  The  material 
properties  for  the  given  problem  are  taken  from  Jackson,  Wagner  and  Wasilewski  [1972] 
and  Qidwai  and  Lagoudas  [2000a]  for  Model  1 and  Set  1 , Brinson  and  Lammering  [1993] 
and  Brinson  [1993]  for  Model  2 and  Set  2,  Ford  and  White  [1996]  for  Model  3 and  Set 
3,  Lexcellent  and  Bourbon  [1996]  for  Model  4 and  Set  4,  Panoskaltsis,  Bahuguna  and 
Soldatos  [2004]  for  Model  5 and  Set  5,  and  Auricchio  [ 1 995]  for  Model  6 and  Set  6.  In  or- 
der to  have  a complete  stress-induced  phase  transformation  from  austenite  and  martensite 
phase,  we  keep  the  temperature  constant  and  applied  the  stress  in  a loading-unloading  cy- 
cle. In  particular  we  consider  three  different  temperature  ranges:  As  < T < Af,T  — Af, 
and  Af  < T undergoing  a cycle  of  loading-unloading.  It  is  recalled  from  the  previous 

17  Since  each  constitutive  model  has  a little  different  approach  for  modeling  and  number  of  parameter 
inputs,  numerical  values  in  parenthesis  are  used  in  the  present  model  to  correlate  those  models. 
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chapter  about  SMA  background  that  at  temperature  T > Af  in  the  stress-free  state,  the 
behaviors  of  SMAs  is  superelastic.  Superelasticity  is  associated  with  ongoing  forward 
and  reverse  phase  transformations  by  a cycle  of  loading-unloading,  in  which  the  observed 
hysteris  effects  offer  evidence  that  the  processes  are  energy-dissipative.  In  the  tempera- 
ture range  T < Aj,  certain  amount  of  permanent  strain  is  observed  after  unloading.  This 
phenomenon  is  due  to  incomplete  reverse  transformation. 

Table  5. 1 . Description  of  material  parameters. 


Material  Description 

Material  Parameters 

Young  modulus 

EA,EM{x  103)  [MPa] 

Poisson  ratio 

uA,  uM 

Thermal  expansion  coefficient 

aA,aM{x KT6)  [/°C] 

Slope  in  stress-temperature  phase  diagram 

CA,CM  [MPa/°C] 

Maximum  transformation  strain 

Martensite  finite  temperature 

Mf°C 

Martensite  start  temperature 

Ms[°  C 

Austenite  start  temperature 

AS[°C 

Austenite  finite  temperature 

u' 

o 

Exponent  in  exponential  hardening  function 

PM  and  (3A 

Thermal  stress  coefficient 

0 [MPa  / °C  ] 

5.1  Isothermal  Uniaxial  Tensile  Loading 

First,  to  characterize  the  evolution  of  the  transformation  strain  and  stress  state  in 
a given  load  increment,  the  simple  uniaxial  tensile  case  is  studied.  Two  elements  are 
used  to  model  along  the  longitudinal  direction  of  a bar.  Because  of  simple  geometry  and 
simple  loading  we  do  not  need  to  use  any  of  EAS  and  ANS  parameters  to  treat  locking 
difficulties  in  usual  FE  formulation.  Figure  5.1  shows  the  geometry  and  applied  loading 
mode  subjected  to  the  uniaxial  tension  considered  here.  The  maximum  applied  force  F 
is  applied  such  that  all  austenite  phase  transformed  into  martensite  phase  producing  a 
maximum  transformation  strain  followed  by  the  full  unloading  to  austenite  state  so  that 
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Table  5.2.  Different  set  of  material  parameters  (Model  1 - Model  3) 


Parameters 

Model  1 (Lagoudas) 

Model  2 (Brinson) 

Model  3 (Ford) 

EA,EM 

70, 30(50) 

67, 23(45) 

53,11(32) 

vA,  uM 

0.30 

0.30 

0.30 

aA,  aM 

11,6.6(8.8) 

8,2(5) 

(5) 

CA,CM 

7.0 

13.8,8.0(11) 

10.0 

<© 

0.05 

0.067 

0.08 

Mf 

-2 

9 

0 

Ms 

18 

18.4 

17 

As 

22 

34.5 

27 

Af 

42 

49 

38 

s 

II 

© 

0.55 

Table  5.3.  Different  set  of  material  parameters  (Model  4 - Model  6) 


Parameters 

Model  4 (Lexcellent) 

Model  5 (Panoskaltsis) 

Model  6 (Auricchio) 

Ea,  Em 

(70) 

67, 23(45) 

1,1 

vA,vM 

(0.30) 

(0.30) 

(0.30) 

aA , aM 

(8.8) 

(8.8) 

(8.8) 

CA , CM 

(6.0) 

(6.5) 

1.0 

£l 

0.064 

(0.05) 

0.10 

Mf 

(-38) 

-33 

10 

Ms 

-23 

-15 

70 

As 

32 

25 

90 

Af 

(47) 

52 

130 

deformed  shape  can  be  recovered  to  the  undeformed  state. 

Material  Set  1. 

Figure  5.3  shows  that  the  input  loading  mode  along  the  load  step,  displacement 
versus  applied  force,  the  evolution  martensite  volume  fraction  along  the  load  step,  the 
stress-strain  response,  the  iteration  number  in  the  global  Newton-Raphson  scheme,  and 
the  transformation  strain  along  the  load  step  at  constant  temperatures  T = 42  = Af 
for  Set  1 . That  is,  we  should  see  the  full  recovery  of  deformed  shape  after  a relatively 
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3-0  structure 


Figure  5.1.  Uniaxial  loading:  Geometry  and  loading. 


barlt.SI20x02x02e1.T42,  square=present,  dashed=  Set  1 (Lagoudas) 


large  deformation,  i.e.,  showing  superelasticity  in  SMA.  As  we  can  see  in  Figure  5.3,  the 
total  number  of  iterations  for  equilibrium  of  force  are  less  than  three  for  both  forward 
and  reverse  transformations,  showing  the  numerical  efficiency  of  the  present  model.  For 
the  case  of  Voce’s  type  exponential  nonlinear  transformation  hardening  of  f3M  = 1.3 
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Uniaxial,  associative  flow  rule,  betaM=betaA=000 
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Figure  5.3.  Uniaxial  loading  (Set  1):  /3M  = (3A  = 0 and  T = 42. 


and  (3A  = 1.5  in  Figure  5.4,  the  total  number  of  iterations  in  Newton-Raphson  scheme 
are  about  six.  Because  of  nonlinearity  of  hardening  function,  iteration  number  in  NR 
schemes  increases.  This  exponential  hardening  curve  in  stress-strain  relation  is  qualita- 
tively similar  to  the  analytic  model  and  experimental  model  for  CuZnAl  SMA  in  Gillet, 
Patoor  and  Berveiller  [1995],  Figure  5.6  shows  the  superelasticity  at  higher  temperature 
than  T > Aj  and  initial  onset  transformation  stress  increases  as  temperature  increases. 
This  plot  is  corresponding  to  process  5 in  Figure  5.5  where  process  5 crosses  over  the 
austenite  start  and  finite  temperature  line  during  reverse  phase  transformation.  Once  the 
specimen  is  fully  recovered,  the  elastic  process  is  carries  over. 

A brief  analytic  procedure  to  identify  the  initial  onset  transformation  stress  and 
saturation  stress  of  transformation  is  now  described.  The  equation  of  line  of  stress  in 
Figure  5.5  reads  as  follows: 

erf  = CM{T-Mf)  and  af  = CM(T  - Ms) 
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Figure  5.4.  Uniaxial  loading  (Set  1):  f3M  = 1.2,  (3 


time  step 

= 1.5,  and  f = 42. 


aA  = Ca(T  - Af)  and  aA  = CA(T  - As)  (5.1) 

so  that  at  T = T (here  T = 42),  dff  = aft1,  aAI  = d^,  aA  = aA,  and  aA  = d^,.  By 
analogy  to  the  well-known  Voce’s  exponential  strain  hardening  (Voce  [1948]),  we  have 

(»"  - So")  = eiC"  (M,  - A/,)  = ~(pAm)  (M,  - M,)  (5.2) 

'I.  (si  - s,(j  = f/.C"'  (.1,  - 4.)  = -t/iA,/,,)  (A;  - .4.)  (5.3) 


where  d^  is  called  the  saturation  stress  at  fully  martensite  transformation  zone  for  £ = 1 
and  Oq1  is  the  initial  onset  martensite  transformation  stress  which  is  corresponding 
to  the  volume  fraction  of  martensite  £ = 0.  Thus,  one  can  calculate  the  initial  onset 
transformation  stress  at  austenite  (or  parent)  state,  i.e.,  when  £ = 0.  For  simplicity  and 
illustration,  let  (3M  = 0 and  m\  = 0,  i.e.,  for  one-dimensional  linear  hardening  model 
case.  The  transformation  function  can  then  be  written 

$(^m)  ;=  aeL  + (pAvo)T-(pA£0)-Y*  = 0 (5.4) 
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Figure  5.5.  Stress-temperature  ( a — T ) phase  diagram  and  stress-strain  curve  of  supere- 
lasticity. 

And  when  £ = 1,  i.e.,  at  fully  martensite  state,  the  transformation  function  is 

aeL  _j_  (pAr^o)  T - (pAf0)  - (~pAr]o)(Ms  - M/)£  - F*  = 0 (5.5) 

Analytical  values  for  ctq7  and  a ^ are  referred  to  the  initial  onset  transformation  stress  at 
168[MPa]  and  the  saturation  stress  (fully  transformed  stress)  at  308 [MPa],  respectively, 
and  can  be  identified  in  Figure  5.3. 

Material  Set  1 : Partial  Recovery. 

Figures  5.7,  5.8,  and  5.9  show  partial  recovery  when  the  specimen  is  tested  at  tem- 
perature T < Af.  The  maximum  applied  force  F is  applied  such  that  all  austenite  trans- 
forms into  martensite  followed  by  the  full  unloading  to  austenite  state.  Figure  5.5  shows 
different  loading  modes  at  different  temperature  in  the  range  of  Ms  < T < Af  and 
Af  < T.  In  each  test,  the  input  loading  mode  along  the  load  step,  displacement  versus 
applied  force,  the  stress-strain  response,  the  evolution  martensite  volume  fraction  along 
the  load  step,  the  iteration  number  in  the  global  Newton-Raphson  scheme,  and  the  stress 
vs.  transformation  strain  are  plotted.  As  we  can  see,  if  the  specimen  is  deformed  at 
above  Ms  but  below  Af  (Ms  < T < A /),  some  stress-induced  martensite  phase,  i.e.,  the 
two  phases  coexist,  will  remain  and  lead  to  an  incomplete  shape  recovery.  In  the  stress- 
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Figure  5.6.  Uniaxial  loading  (Set  1):  = (3A  = 0 and  Af  < T = 50. 


temperature  phase  diagram  in  Figure  5.5,  the  loading  mode  cases  2 and  3 will  represent 
this  partial  recovery  stress-induced  transformation  response. 

Material  Set  2 - Set  6. 

One  of  the  pioneering  constitutive  model  developed  by  Brinson  and  Lammering 
[1993]  has  been  widely  used  because  of  its  simple,  wide  range  of  transformation  temper- 
atures, as  well  as  their  FE  formulation  of  truss  element,  which  seems  to  be  the  first  FE 
formulation  shown  explicitly  for  shape-memory  materials.  They  have  also  compared  with 
experimental  input  data  which  is  used  by  cosine  model  of  Liang  and  Rogers  [1990],  The 
stress-strain  curve  for  superelastic  effect  is  demonstrated  in  Figures  5.10,  5.1 1,  and  5.12 
at  temperatures  T — 40,  T = 50,  and  T = 60.  For  each  test  of  Material  Set  2,  the  present 
model  is  compared  to  the  results  from  Model  2 where  Material  Set  2 is  used. 

In  Ford  and  White  [1996]  model,  Brinson’s  model  is  modified  to  allow  prediction 
of  mechanical  behavior  into  region  for  plasticity.  To  accomplish  this  it  is  assumed  that 
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Figure  5.7.  Uniaxial  loading  (Set  1):  (3M  = (3A  = 0 and  As  < T = 32  < Af. 


after  yielding  SMA  behaves  as  a linear  work-hardening  plastic  material.  The  mechani- 
cal testing  was  done  to  acquire  uniaxial  stress-strain  data  for  SMA  wires  up  to  failure. 
Figures  5.13  and  5.14  showed  that  our  model  is  valid  at  temperature  T & Af. 

Lexcellent  and  Bourbon  [1996]  presented  experimental  data  for  Ni-Ti  and  Cu-Zn- 
A1  polycrtstalline  materials  and  compared  with  the  results  from  their  constitutive  models 
(Figures  5.15  and  5.16).  The  present  model  is  also  be  compared  to  their  model,  as  well  as 
their  experimental  data  for  Ni-Ti  system  of  SMAs. 

Recently,  Lubliner,  Taylor  and  Auricchio  [1993]  and  Auricchio  and  Taylor  [1994] 
developed  generalized  plasticity  to  effectively  describe  the  ratchetting  and  real  material 
response  in  elasto-plasticity.  The  response  of  shape-memory  materials  has  also  been  eas- 
ily described  by  generalized  plasticity  (Auricchio  [1995]).  The  theoretical  framework  by 
Panoskaltsis  et  al.  [2004]  for  SMA  model  is  also  based  on  generalized  plasticity  theory. 
See  figures  5.17,  5.18,  5.19,  and  5.20  for  comparison  results  with  the  present  model.  Note 
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Figure  5.8.  Uniaxial  loading  (Set  1):  f3M  = (3A  = 0 and  As  < T — 30  < Af. 
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Figure  5.9.  Uniaxial  loading  (Set  1):  (3M  = /3A  = 0 and  As  < f = 25  < Af. 
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that  the  Material  Set  6 has  been  referred  to  hypothesis  alloy.  As  we  can  see,  our  model 
seems  to  reproduce  other  models  effectively  and  also  be  compared  with  experimental 
data  at  around  austenite  finite  temperature  T « Af.  This  will  support  the  validation  of 
our  model. 


Figure  5.10.  Uniaxial  loading  (Set  2):  = @A  = 0 and  Af  < T = 60.  Comparison 

between  Model  2 with  dotted  line  and  present  model  with  squared  line. 


5.2  Simply  Supported  Beam 

Consider  a simply-supported  beam  with  a rectangular  cross-section  subjected  to  a 
concentrated  load  at  the  mid-span  of  the  beam  to  simulate  a condition  of  pure  bending, 
as  shown  in  Figure  5.21.  Figure  5.22  shows  FE  modeling  of  10  elements  in  longitudinal 
direction  and  1 element  through  the  thickness  direction  of  beam.  The  purpose  is  to  de- 
termine the  stresses  and  deflections  in  the  mid-span  of  the  beam  in  the  superelastic  state. 
For  the  three-point  bending,  the  plate  is  simply  supported  on  two  edges  and  free  at  the 
other  two  edges  and  the  load  consists  of  one  line  load  at  the  mid-span  of  the  beam.  The 
beam  has  a rectangular  cross  section  of  a width  of  1 .0  and  a thickness  of  1 .0.  The  length 
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barU.SI20x02x02e2.T50,  square=present,  dashed=  Set  2 (Brinson) 


Figure  5.11.  Uniaxial  loading  (Set  2):  (3M  = /3A  = 0 and  Af  < T = 50.  Comparison 
between  Model  2 with  dotted  line  and  present  model  with  squared  line. 


barU.SI20x02x02e2.T40,  square=present,  dashed=  Set  2 (Brinson) 


Figure  5.12.  Uniaxial  loading  (Set  2):  (3M  = (3A  = 0 and  As  < T = 40  < Af. 
Comparison  between  Model  2 with  dotted  line  and  present  model  with  squared  line. 
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barU.SI20x02x02e3.T45,  square=present,  dashed=  Set  3 (Ford) 


Figure  5.13.  Uniaxial  loading  (Set  3):  /3M  = (3A  = 0 and  Af  < T = 45.  Comparison 
between  Model  3 with  dotted  line  and  present  model  with  squared  line. 


barU.SI20x02x02e3.T32,  square=present,  dashed=  Set  3 (Ford) 


Figure  5.14.  Uniaxial  loading  (Set  3):  (3M  = f3A  = 0 and  As  < T — 32  < Af. 
Comparison  between  Model  3 with  dotted  line  and  present  model  with  squared  line. 


104 


barU  e7.100,  square=present,  dashed=set  7,  star  = exp 


Figure  5.15.  Uniaxial  loading  (Set  4):  /3M  = (3A  = 0 and  Af  < f = 100.  Comparison 
between  Model  4 with  dotted  line  and  present  model  with  squared  line. 


barU  e7.80,  square=present,  dashed=set  7,  star  = exp 


Figure  5.16.  Uniaxial  loading  (Set  4):  (3M  = (3A  = 0 and  Af  < T = 80.  Comparison 
between  Model  4 with  dotted  line  and  present  model  with  squared  line. 

between  the  two  supports  is  20.0.  Due  to  symmetry,  only  half  of  the  beam  is  modeled 
in  FE  analysis  where  ten  elements  are  meshed  along  the  longitudinal  direction  of  beam. 
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barU.SI20x02x02e5.24,  square=present,  dashed=set  5 
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Figure  5.17.  Uniaxial  loading  (Set  5):  f3M  = f3A  = 0 and  T = 24  < A/.  Comparison 
between  Model  5 with  dotted  line  and  present  model  with  squared  line. 


The  material  properties  are  the  same  as  those  used  in  uniaxial  tensile  case  (Table  5.2). 
Figure  5.23,  5.24,  and  5.25  show  the  displacement  response  at  mid-center  where  the  load 
is  applied.  These  plots  show  the  numerical  efficiency  of  the  proposed  SMA  constitutive 
model  where  the  number  of  global  Newton-Raphson  iteration  steps  along  the  load  step  is 
six  that  is  obtained  uniformly  along  the  whole  load  increments.  Figure  5.26  shows  de- 
formed shape.  In  order  to  see  the  stress  distribution  through  the  thickness  of  the  beam, 
we  may  need  to  stack  a lot  of  element  in  the  thickness  direction  as  shown  in  Figure  5.27 
where  10  elements  in  longitudinal  direction  and  12  elements  through  the  thickness  direc- 
tion of  beam  are  discretized.  A beam  of  rectangular  cross-section  of  a width  of  1 .0  and  a 
thickness  of  1 .0,  the  length  of  4.0  between  the  two  supports  is  considered. 

In  Figures  5.28  and  5.29,  the  normal  stress  distribution  are  plotted  in  a three-point 
bending  test.  Because  of  volume  conserving  effects  during  phase  transformation  simi- 
lar to  plasticity,  incompressibility  effects  is  preserved  by  using  enhanced  assumed  strain 
mode  where  the  number  of  EAS  parameter  is  seven.  Due  to  symmetry,  only  half  of  the 
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barU.SI20x02x02e5.36,  square=present,  dashed=set  5 


Figure  5.18.  Uniaxial  loading  (Set  5):  /3M  — (3A  = 0 and  T = 36  < Af 
between  Model  5 with  dotted  line  and  present  model  with  squared  line. 


barU.SI20x02x02e5.48,  square=present,  dashed=set  5 


Figure  5.19.  Uniaxial  loading  (Set  5):  (3M  = (3A  = 0 and  T = 48  < Af. 
between  Model  5 with  dotted  line  and  present  model  with  squared  line. 
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barU  e6.160,  square=present,  dashed=set  6 
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Figure  5.20.  Uniaxial  loading  (Set  6):  /3M  = /3A  = 0 and  Af  < T = 160.  Comparison 
between  Model  6 with  dotted  line  and  present  model  with  squared  line. 
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Figure  5.21.  Three-point  bending  problem:  Geometry  and  loading. 

beam  is  modeled  in  FE  analysis.  During  the  forward  phase  transformation,  the  stress 
distribution  in  the  cross-section  evolves  as  follows: 

• Increasing  the  intensity  of  the  load,  the  beam  will  be  in  a fully  elastic  state  as  long  as 
the  maximum  normal  stress  does  not  reach  the  initial  onset  of  phase  transformation 
stress.  That  is,  no  phase  transformation  occurs  in  the  cross-section.  Thus,  the  range 
of  martensite  volume  fraction  is  f < 0. 
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3-D  structure 


Figure  5.22.  FE  model  of  three-point  bending  beam:  1 element  through  thickness. 

• As  the  load  increases  above  the  elastic  limit  load,  phase  transformed  zones  de- 
velop in  the  upper  and  lower  portions  of  the  beam.  Like  in  the  plasticity,  the  initial 
transform  occurs  at  the  outermost  fibres  on  the  central  section  when  the  maximum 
applied  moment  reaches  the  moment,  corresponding  to  the  stress  at  the  onset  of 
transformation  in  loading,  for  the  the  beam  cross-section.  Two  symmetric  marten- 
sitic transformed  zones  spread  inward.  In  these  regions,  the  range  of  martensite 
volume  fraction  is  0 < £ < 1. 

• Distinguished  material  response  that  is  quite  different  from  plasticity  occurs  when 
£ > 1,  i.e.,  when  the  phase  transformation  is  completed.  We  can  see  a material 
re-stiffening  in  that  area  is  encountered  in  that  regions  where  the  material  behaves 
like  an  elastic  response  again.  The  range  of  martensite  volume  fraction  is  £ = 1 

Figure  5.30  shows  the  same  responsemf  three-point  bending  problem  when  we  use 
pure  displacement-based  20-node  solid  element.  In  the  case  of  a concentrated  force,  these 
boundary  lines  are  parabolas  while  for  the  case  of  a uniformly  distributed  load,  we  would 
have  to  see  that  the  transformed  rezones  are  bounded  by  straight  lines,  like  in  the  case  of 
plasticity. 
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Figure  5.23.  Three-point  bending:  Based  on  8-node  solid-shell,  ANS=2,  EAS=4. 
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Figure  5.24.  Three-point  bending:  Based-on  8-node  solid-shell,  ANS=0,  EAS=0 
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Figure  5.25.  Three-point  bending:  Based-on  20-node  displacement-based  solid  element. 
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Figure  5.26.  Three-point  bending:  Undeformed  and  deformed  shape  of  ANS(2)/EAS(2) 
turn  on/off. 
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3-D  structure 


Figure  5.27.  FE  model  of  three-point  bending  beam:  12  elements  through  thickness. 


sp3c.SI02x01  xl  eM0.2b1 0x1 2L2.fi  40.202.Sxx.G02to50 


Figure  5.28.  Three-point  bending  problem.  Normal  stress  (axx)  distribution  at  gaussian 
integration  point.  Based-on  8-node  solid-shell  element. 
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sp3c.SI02x01  xl  eM0.2b1 0x1 2L2.fi  40.202.Sxx.N02to50 


Figure  5.29.  Three-point  bending  problem.  Normal  stress  ( axx ) distribution  at  nodes. 
Based-on  8-node  solid-shell  element. 


sp3c.SI02x01  xl  eM0.04x08.f140.Sxx:  G02to50 


Figure  5.30.  Three-point  bending  problem.  Normal  stress  ( axx ) distribution  at  gaussian 
integration  point.  Based-on  20-node  displacement-based  solid  element. 
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5.3  Asymmetry  in  Tension  and  Compression 


The  non-symmetry  between  tensile  and  compressive  tests  experimentally  observed 
in  shape-memory  alloys  is  taken  into  account  through  the  definition  of  non-symmetric 
transformation  function  (Orgeas  and  Favier  [1995]).  In  order  to  characterize  the  asym- 
metry effect  subject  to  the  tension  and  compression  during  the  phase  transformation,  a 
simple  uniaxial  loading-unloading  cycle  is  applied  to  the  specimen.  Only  one  element 
is  used  to  model  with  8-node  solid-shell  element  and  it  is  not  necessary  to  turn  on  the 
EAS  and  ANS  parameters  because  of  simple  geometry  and  simple  loading  mode  and  the 
FE  model  is  the  same  as  in  Figure  5.1.  The  loading-unloading  cycle  to  test  the  asym- 
metry effect  is  shown  in  Figure  5.31.  The  maximum  applied  force  F is  applied  such 
that  all  of  the  austenite  transforms  into  martensite,  producing  a transformation  strain  of 
5%  followed  by  the  full  unloading  to  austenite  state  for  regular  symmetry  in  tension  and 
compression.  Then  the  specimen  is  compressed  in  forward  transformation  and  released 
to  original  shape  in  reverse  phase  transformation  (Figure  5.31).  Recall  the  transformation 
function  of  J2  — J3  flow  from  the  previous  section: 


where  ac  is  the  initial  yield  stress  in  compression  and  at  is  the  initial  yield  stress  in 
tension.  When  b — 0,  the  function  will  reduced  to  von-Mises  type  J2  flow  yield  function. 
If  we  let  ac  = k x at,  then 


(5.6) 


where 


(5.7) 


Vrf  k2  - 1 


(5.8) 


so  that  k = 1 will  be  reduced  to  b — 0 and  k — 1.40  will  yields  b = 0.85  (Gillet  et  al. 
[1998],  Terriault  et  al.  [1999]). 
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Uniaxial,  associative  flow,  Asymmetry,  b=1 .0,  betaM=betaA=0 


time  step  time  step 


Figure  5.31.  Symmetry  in  uniaxial  loading:  k = 1.00  and  (3M  = /3A  = 0. 

Figures  5.31, 5.32,  and  5.33  read  the  input  applied  loading  mode  along  the  time,  the 
stress-strain  response,  the  global  Newton-Raphson  iteration  numbers,  and  the  total  strain 
along  the  time  at  temperatures  of  T = 42,  i.e.,  at  austenite  finite  temperature.  As  we  can 
see  in  Figure  5.32,  the  total  number  of  iterations  being  satisfied  the  equilibrium  of  force 
is  less  than  three  for  both  forward  and  reverse  transformations.  For  the  case  of  nonlin- 
ear transformation  hardening  of  (3M  — 1.3  and  (3A  = 1.5,  Figure  5.33  shows  the  total 
number  of  iterations  in  Newton-Raphson  scheme.  Although  experimental  results  showed 
that  there  is  a significant  asymmetry  between  and  compression  during  deformation  due 
to  martensitic  phase  transformation,  the  precise  origin  of  the  asymmetry  is  not  known  yet 
and  is  to  be  more  closely  investigated  for  its  micro-structural  origins.  The  origin  of  such 
an  asymmetry  is  neither  due  to  the  anisotropy  of  material,  nor  to  a Bauschinger  effect  in- 
duced by  mechanical  cycling  load,  nor  to  the  material  being  sensitive  to  hydrostatic  pres- 
sure.18 This  asymmetry  may  be  due  to  different  modes  of  deformation  accommodation 

18  Auricchio  et  al.  [1997]  has  used  I\  — J2  Drucker-Prager  yield  function  to  accommodate  the  volume 
change  during  phase  transformation  and  showed  asymmetry  response  in  tension  and  compression. 
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Uniaxial,  associative  flow,  Asymmetry,  b=1.05,  betaM=betaA=0 


time  step  time  step 


Figure  5.32.  Asymmetry  in  uniaxial  loading:  k = 1.05  and  (3M  = (3A  = 0. 

during  stress-induced  martensite  phase  transformation  (SIMPT).  Thus,  from  a mechani- 
cal standpoint,  it  appears  that  a correct  modeling  of  the  Ni-Ti  response  has  to  take  into 
account  the  role  of  the  third  invariant  of  the  stress  tensors,  as  it  has  already  been  proposed 
by  Gillet  et  al.  [1998]  on  a Cu-Zn-Al  shape-memory  alloy,  using  a Prager  criterion  to 
characterize  phenomenologically  the  tension-compression  asymmetry. 

5.4  Pinched  Cylindrical  Shell 

This  problem  was  widely  tested  in  elasto-plastic  problem  as  one  of  benchmark  test 
of  shell  element  formulation  (Vu-Quoc  and  Tan  [2003];  Li,  Hao  and  Liu  [2000]).  A 
cylinder  is  loaded  by  an  inward  radial  load  acting  at  the  middle  section  of  the  cylinder 
(Figure  5.34).  At  both  ends  the  boundary  conditions  have  been  fixed  such  that  the  circular 
shape  of  the  end  cross  section  is  preserved  and  free  deformation  in  axial  direction  is 
allowed.  The  material  properties  are  the  same  as  the  Material  Set  1 in  Table  5.2  and 
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Uniaxial,  associative  flow,  Asymmetry,  b=1.05,  betaM=1.3,  betaA=1.5 


time  step  time  step 


Figure  5.33.  Asymmetry  in  uniaxial  loading:  k = 1.05,  j3M  = 1.3,  and  (3A  = 1.5. 
geometric  properties  are 

R = 300,  H = 600,  h = 0.03  (5.9) 

where  R is  the  radius,  H the  height,  and  h the  thickness  of  the  cylinder. 

Because  of  the  symmetry,  only  one-eighth  of  the  cylinder  needs  to  be  modeled  by 
16  x 16  x 1 solid-shell  elements  with  578  nodes  (Figure  5.34).  For  the  given  problem, 
maximum  magnitude  of  load  f=300000  is  incrementally  applied  in  50  load  steps  to  finish 
the  run.  The  deformed  shape  of  cylinder  at  load  step  10,  20,  40,  and  50  are  shown  in 
Figures  5.35  and  5.36.  As  one  of  reference,  we  test  elastic  material  with  E = 50  x 103 
and  v — 0.3  and  same  number  of  load  step  divided  into  50  shown  in  Figure  5.37.  At  the 
end  of  our  computation,  the  two  opposite  points  of  inner  surface  of  the  cylindrical  shell 
come  together  appear  more  buckling  modes  than  elastic  material.  It  is  observed  that  the 
response  from  the  linear  elastic  material  and  the  shape-memory  material  are  very  close, 
which  justify  the  use  of  the  present  solid-shell  element  for  the  shell  structures. 
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Figure  5.34.  Large  motion  of  pinched  cylindrical  shell:  Geometry  and  loading. 


Pcyl.SI.256e3.f30000.204.t01 0 


Pcyl.SI.256e3.f30000.204.t020 


Figure  5.35.  Pinched  cylinder:  Deformation  at  load  step  10  (left)  and  load  step  30  (right). 
ANS=2,  EAS=4. 
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Pcyl.EL256e3.f30000.204.t040 


Pcyl.SI.256e3.f30000.204.t050 


Figure  5.36.  Pinched  cylinder:  Deformation  at  load  step  40  (left)  and  load  step  50  (right). 
ANS=2,  EAS=4. 


Pcyl.EL256e3.f30000.204.t030 


Pcyl.EL256e3.f30000.204.t050 


Figure  5.37.  Pinched  cylinder:  Deformation  at  load  step  30  (left)  and  load  step  50  (right) 
ANS=2,  EAS=4,  Elastic  material. 


CHAPTER  6 
CONCLUSION 


This  study  addresses  a simple  but  numerically  efficient  three-dimensional  unified 
thermodynamic  constitutive  model  for  the  analysis  of  superelastic  response  of  shape- 
memory  alloys.  To  assess  the  viability  of  the  proposed  constitutive  model,  some  simple 
and  typical  SMA  stmctures  such  as  uniaxial  test,  three-point  bending  as  well  as  asymme- 
try in  tension  and  compression  were  simulated.  The  implementation  of  the  SMA  model 
in  a solid-shell  finite  element  framework  is  also  addressed.  In  order  to  remedy  locking 
problems  encountered  in  low-order  finite-element  formulation,  we  use  solid-shell  element 
with  an  optimal  choice  of  EAS  parameters. 

Numerical  examples  confirm  that  the  proposed  constitutive  model  performed  much 
better  than  the  unified  thermodynamic  model  in  the  sense  that  the  number  of  global 
Newton-Raphson  iteration  scheme  can  be  reduced  significantly.  The  proposed  model 
is  also  comparable  with  classical  plasticity  in  the  sense  of  being  an  additive  decompo- 
sition of  free  energy  and  associative  flow  rule  and  normality  in  time  continuous  model 
and  closes-point  projection  method  in  algorithmic  (time  discrete)  model.  This  numerical 
performance  leads  to  the  conclusion  that  the  present  method  is  an  effective  and  simple 
computational  tool  for  analysis  of  superelastic  SMA-based  devices. 

As  we  can  see  in  the  previous  section,  our  constitutive  model,  which  is  based  on 
free  energy  of  function  of  strain  and  temperature  with  which  most  of  plasticity  constitutive 
model  has  been  developed,  can  easily  be  extended  to  describe  a conventional  plasticity 
zone  which  is  beyond  the  transformation  zone  in  shape-memory  alloys. 
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APPENDIX  A 

DERIVATIONS  OF  SOME  USEFUL  FORMULAE 


A.  1 Derivation  of  Gradient  of  Transformation  Direction  Tensor 

We  show  detailed  mathematical  procedures  for  derivation  of  gradient  of  transfor- 
mation direction  tensor  A in  case  of  J2  flow  used  in  implicit  backward  Euler  rule  in  stress 
integration  scheme. 


(A.l) 


where  the  transformation  tensor  A is  assumed  to  have  the  following  form 

a = { <i>°> 

\ O/^pT  € r > (I  < 0) 

and  the  equivalent  stress  o'  is  defined  as 

1/2 


(A. 2) 


(A. 3) 


Here  I is  the  4th-order  symmetric  identity  tensor  and  its  indicial  notation  is  given  by 


I ^ d-  f l ^ j ® 


(A.4) 


and  1 is  the  2nd-order  identity  tensor  given  by 


1 — Sij  6j  A'1  Gj 


(A. 5) 


We  begin  with  computing  the  partial  derivative  of  deviatoric  stress  with  respect  to 
the  stress  as  follows: 


da'ij  doij  1 dakk  1 . . . 1 
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(A. 6) 
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120 


121 


d 

A,  a' 

1 

1 3 d{a'ija'ij) 

dakt 

y 2 y 

2> 

l^o'  a'  2 dcr^ 

/ 2 PQ  PQ  Kl 

1 1 3 

2 ./-( j'  (j'  2 

V 2 upqupq 


daij  / , / d<i 

—Va  + <?„ L 


da , 


kl 


y 


y 


1 1 3 ro  , . 3 o£, 

2 TT7 — 7~  2 ~ 9 /3  / / 

2 V 2^9  2 2 V^pX 


31  , 

2a'  kl 


PQ 


Defining  5 :=  \a'pqa'pq,  the  above  equation  can  be  rewritten  as 


2 dA 
3eL  da 


1 da\,  d 


s/Bdc 


pq 


13  da. 


pq 


i , , s(3/24«tm-'/2 


\Z~B  da 


pq 


i 


+ %• 


dcr 


PQ 


y/B  da 


pq 


or'  ^3  , , \—3/2 

aij  2 ( 2 ^b^y  ' 


^.+ 

<9an 


/ d<j 


'pq 


1 da' 


y 


1 


1 


y/B  da. 


a 


y 


pq 


1 5(Ty 


2 *-2 


1(20 


, 3 


cr,„- 


y/Bdapq  ^2  3 , , /3^— 7 

^ 2 yuy  V 2uvu*. 

[2  1 


=cr. 


P9 


y 


<9  A 

3 cr 


'5<Ty 

_/  _/  ' 
7ijapq 

/ 

ij 

d&pq 

/! 

1 

<jaM 

/ 3 

dapq 

a'ija'ij 

1 da, 


PQ  / 


1 

dapq 

- / 

X. 

X. 

1,  X 
2 °ij  °pq 

<3a'pQ 

W 

vipvjq 

1 cr ' 

3 cr ' 


The  above  expression  is  used  in  calculating  of  algorithmic  tangent  modulus. 


(A. 8) 


(A. 9) 


(A.  10) 
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A. 2 Derivation  of  Incremental  Stress 


For  continuum  tangent  moduli  (4.107),  some  tensor  algebra  calculation  is  neces- 
sary to  determine  the  incremental  stress  d a . Here  we  show  detailed  mathematical  pro- 
cedures:19 


d a 


...  + 


-J  ^ \ 1 
da  J a 


a dT 


(C  : r) 
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: e 


df_ 
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For  the  right-handed  side  terms  being  in  an  indicial  notation, 
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19  Most  of  the  derivation  is  very  similar  to  the  case  of  plasticity.  See  the  lecture  note  of  Nonlinear  Finite 
Element  Method,  EGM6934,  1998. 
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For  the  left-handed  side  terms  being  in  an  indicial  notation, 
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By  comparing  (A.  17)  with  (A.  19),  the  proof  has  been  shown. 


(A.  18) 


(A.  19) 


A.3  Recover  1-D  quantity  from  3-D  of  Transformation  Direction  Tensor 


For  the  case  of  J2  flow,  transformation  direction  tensor  A provides  the  directions 
in  which  the  transformation  strains  develop.  For  example,  the  transformation  strain  is 
going  to  develop  along  the  direction  of  the  stress  during  forward  transformation  (£  > 0), 
while  it  is  going  to  develop  along  the  opposite  direction  of  the  previously  developed 
transformation  strain  during  reverse  transformation  (£  < 0).  In  this  section,  we  show  that 
how  the  transformation  direction  tensor  A from  3-D  to  1-D  quantity. 


124 


Equivalent  stress,  (a). 


ItLjS,  (£  > 0) 
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(A. 20) 


eL£e?,  (£  < 0) 

cr 

During  the  reverse  transformation,  all  of  the  variants  must  return  to  the  same  parent  sym- 
metry. Thus,  the  reverse  form  of  A should  be  governed  by  kinematical,  or  symmetry 
considerations  (Boyd  and  Lagoudas  [1996,  p.820]). 
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Assume  an  = ay  for  uniaxial  loading  in  1 -dimensional  case. 

\ay  0 

0 
0 


&ij 


ay  3 ay 


i 

' 3 ay 


’y 


2 '2 
3°y . 


0 

r 2 
3 °y 

0 

0 

0 

= 

0 

~\°y 

0 

. 

0 

0 

l _ 

3 ay  . 

/ 1 

\2 

1 \2 

2 

+ (~3 

<7y) 

+ (- 

3 V 

= 3W 

a = 


-s  : s 


(A.21) 


(A.22) 


(A. 23) 


— | 2 11  ~ a22)2  + (&22  — 033  ) 2 + (^33  ~ 0"ll)2  +3  ^23  + <t|j  + a\2  | 


1/2 


= V§§M2  = *» 


Thus  for  1 -dimensional  case 
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Equivalent  transformation  strain,  (ept).  Recall  equivalent  transformation  strain  similar 
to  plastic  strain  in  plasticity  theory: 
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In  hardening  hypothesis,  it  is  assumed  that  the  extent  of  hardening  is  a function  of  the 
total  equivalent  transformation  strain.  Considering  volume  conservation  during  inelastic 
flow,  we  obtain 
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(A. 26) 


e*  = 


(rf  )2  + (ef)2  + (ef  )2 
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(A. 27) 


Note  that  if  the  SMA  wire  is  loaded  under  tension  only  for  1-d  case,  we  can  reduce 
the  transformation  direction  tensor  A is  a positive  constant  equal  to  eL. 


APPENDIX  B 

MAXIMUM  TRANSFORMATION  DISSIPATION  PRINCIPLE. 


The  principle  of  maximum  inelasticity,  e.g.,  plasticity  or  phase  transformation  in 
solids,  plays  a crucial  role  in  the  variational  formulation  of  inelasticity  problems  such 
as  plasticity  and  phase  transformation  in  solids  (Simo  and  Hughes  [1998,  p.98]).  The 
fundamental  significance  of  this  principle  implies: 


• Associative  flow  rule  of  the  elastic  domain  in  stress  space, 

• Normality  of  flow  rule  of  the  elastic  domain  in  stress  space. 


• Convexity  of  the  elastic  domain  C. 


The  mathematical  theory  of  convex  analysis  provides  an  elegant  framework  to  pose 
the  alternative  forms  of  the  flow  law.  Here  it  is  now  necessary  to  introduce  the  existence  of 
an  elastic  region  and  transformation  surface,  which  then  goes  on  to  deduce  the  convexity 
of  the  elastic  region  and  normality  of  the  inelastic  strain  rate  to  the  transformation  surface 
(Reddy  and  Martin  [1994,  p.435]).20  For  example,  in  the  classical  theory  of  plasticity, 
convexity  of  the  elastic  region  and  the  normality  law  are  consequences  of  the  maximum 
plastic  work  inequality. 

Let  us  recall  the  intrinsic  dissipation  inequality  that  we  used  in  the  previous  section 
(4.24): 


font  ■=  V ■ ept  - > o 


(B.l) 


We  can  rewrite  the  above  equation  in  the  form  of  a scalar  product  as  the  function  </>(X) 
which  is  differentiable  for  X f 0 and  is  expressed  as  follows: 


<f>(X)  = Y-X>0 


(B.2) 


20  For  a function  / from  A to  B such  that  / : A — > B,  we  call  set  A the  domain  of  the  function  and  set 
B the  range  (or  region)  of  the  function. 
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where 

dJ- 

Y = {cr,-p—)  and  X=(e^|)  (B.3) 

called  the  generalized  stresses  and  the  generalized  strain  rates,  respectively.  The  general- 
ized stresses  Y are  not  unbounded  and  they  are  restrained  inside  a convex  elastic  domain 
C which  contains  the  origin  O (Maugin  [1992,  p.44]). 

Remark  B.l.  In  the  classical  theory  of  an  elastic-perfectly  plastic  material  for  which 
case  Y = cr  (Simo  and  Hughes  [1998,  p.98]),  the  principle  of  maximum  plastic  dissipa- 
tion states  that,  for  given  plastic  strain  rate  ep  among  all  possible  stresses  r satisfying 
the  yield  function  (i.e.,  fy(r)  = 0),  the  plastic  dissipation  (or  the  rate  of  plastic  work) 
which  is  now  given  by 


<j)(  ep)  = cr  : ep 

attains  its  maximum  for  the  actual  stress  tensor  cr  , that  is, 


(B.4) 


cr  : cp>  t : ep 


(B.5) 


for  all  admissible  stresses  r,  i.e.,  stresses  satisfying  fv(T)  < 0-  Alternatively,  the 
postulate  of  maximum  plastic  work  states  that 


< Kate(ep ) = max 
t eC 


/»(t)  < 0 


(B.6) 


that  is  the  maximum  being  taken  over  all  r such  that  fy(r)  < 0. 

With  the  additional  demand  of  associated  plasticity  that  the  direction  of  the  plastic 
strain  rate  must  be  orthogonal  to  the  yield  surface,  the  principle  of  maximum  dissipation 
states  that  the  actual  stress  state  cr  maximizes  the  dissipation  for  all  possible  stress  states 
r 6 C: 


= a 


ep  — 


f*{cr)~ 


t : ep  > 0 V r G C. 
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Figure  B.2  shows  a geometric  interpretation  of  this  principle.  Choosing  r = 0 shows 
that  the  dissipation  must  be  non-negative  for  all  cr  £ C such  that  cr  : e p > 0. 

The  principle  of  maximum  dissipation  can  now  in  turn  be  used  for  the  construction 
of  a normal  rule,  describing  the  evolution  of  the  plastic  flow.  For  a given  convex  domain 
the  stress  space  C,  the  associated  flow  rule  can  be  found  as  the  solution  of  an  optimization 
problem  which  can  be  written  as  a saddle-point  problem  using  the  appropriate  Lagrange 
function  (Simo  and  Taylor  [1985,  p.98]) 


where  7 is  called  the  plastic  multiplier.  From  this  and  the  Kuhn-Tucker  conditions 


£(  o- , 7)  :=  7 /„(  cr,  q) -<!>*(  or) 


(B.8) 


7 > 0,  fy  < 0,  7 fy  = 0 


(B.9) 


we  obtain  the  flow  rule 


(B.  10) 


Likewise,  the  evolution  rule  for  hardening  can  also  be  obtained 


(B.ll) 


where  a represents  a hardening  variable. 


I 


Mathematically  speaking,  dissipation  inequality  is  expressed  then  in  terms  of  X in 


the  form 


(Y*ec) 


(B.12) 


The  dissipation  inequality  (B.12)  is  equivalent  to  the  inequality 


Y-X  > Y* *X 


(B.13) 
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holds  for  all  admissible  generalized  stresses  Y*  G C.  In  another  form,  the  above  equation 
can  be  stated  as  (Figure  B.2) 


(Y  — Y*) *X  > 0,  V Y*  G C 


(B.14) 


which  is  called  the  Hill-Mandel  maximal-dissipation  principle  (Hill  [ 1 948],  Maugin  [ 1 992, 
p.l  15]).  In  particular,  we  see  from  (B.l)  that  (B.14)  is  in  accord  with  the  reduced  dissipa- 
tion inequality  since  we  have,  choosing  Y*  = 0 (remember  that  0 G C by  assumption). 


Y*X  > 0 <=>  tr  : ept+  A"  : otv  > 0. 


(B.15) 


Elastic  Domain. 

The  existence  of  a region  TZ  in  the  space  of  generalized  stresses,  called  the  elastic 
region,  is  assumed.  This  region  has  the  property  that  elastic  behavior  occurs  at  a material 
point  at  which  Y G TZ.  The  boundary  of  this  region  is  denoted  by  dlZ,  and  is  referred 
to  as  the  transformation  surface.  It  is  often  convenient  to  represent  the  transformation 
surface  by  the  level  surface  (Figure  B.4)  of  a function  f*,  so  that 


dlZ  = 


{y|^.*(Y)  = o} 


(B.l  6) 


where  we  have  chosen  f*  in  such  a way  that  its  level  set  at  zero  coincides  with  <97 Z.  The 
elastic  region  is  the  set 


7Z  = 


{y|^*(Y)  < o} 


(B.17) 


The  union  of  the  elastic  region  with  the  transformation  (yield)  surface  together  make  up 
the  region  of  admissible  generalized  stresses,  the  so-called  elastic  domain,  denoted  by  C 


C = TZU  dlZ  or  C = { Y 


<t>*( Y)  < 0 


(B.l  8) 


so  that  the  transformation  (yield)  surface  is  the  set  of  points  Y that  satisfy  </>*(Y)  = 0, 
while  the  elastic  region  is  given  by  the  set  of  points  Y for  which  0*(Y)  < 0 (Figure  B.l). 
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Figure  B.l.  Transformation  (or  yield)  surface  and  elastic  region. 


Generalized  Normality. 

There  are  two  major  consequences  of  the  maximum  inelastic  work  inequality.  First, 
it  can  be  shown  that  the  principle  implies  the  existence  of  a closed  convex  set  C of  achiev- 
able values  of  the  generalized  stress,  so-called  elastic  domain  C.  Secondly,  it  can  be 
shown  that  the  generalized  inelastic  strain  rate  X associated  with  a generalized  stress  Y 
on  the  transformation  (yield)  surface  dTZ  is  normal  to  the  tangent  hyperplane  at  the  point 
Y to  the  transformation  (yield)  surface  dTZ.  This  result  is  generally  referred  to  as  the 
generalized  normality  rule.  In  the  event  that  the  transformation  (yield)  surface  is  non- 
smooth, the  normality  law  in  the  framework  of  Generalized  Standard  Materials  (GSM) 
states  that  the  generalized  inelastic  strain  rate  X is  a vector  belonging  to  the  cone  of  the 
outward  normals  in  Y at  the  elastic  domain  C, 


X G Nc( Y) 


(B.19) 


which  may  be  considered  as  an  evolution  law  for  the  generalized  inelastic  strain  X.21  (see 
an  illustration  in  Figure  B.2).  Nc{ Y)  represents  the  cone  of  the  outward  normals  in  Y at 
the  elastic  domain  C.  The  above  condition  corresponds  to  the  fact  that  the  angle  between 

21  For  the  case  of  plasticity,  the  maximal-dissipation  principle  (B.14)  is  expressing  the  following  nor- 
mality property:  for  all  a * to  the  inward  side  of  the  tangent  to  the  yield  surface  at  er , then  e p must  be 
directed  along  the  outward  normal  there. 
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Figure  B.2.  Normal  cone  to  a convex  set:  Smooth  surface  (a)  and  non-smooth  surface 
case  (b). 


YY*  and  X is  obtuse  (>  7r/2).  The  normal  cone  to  C is  also  defined  by 


NC(Y)  :=  jx  [Y  — Y*]*X  > 0 , VY*ecj 


(B.20) 


Flow  Rule  and  Consistency  Condition. 

If  Y is  interior  to  C,  then  the  response  is  purely  reversible,  i.e.,  purely  elastic.  That 
is, 


Y G int(C)  =»  Nc( Y)  = {0} 


(B.21) 


where  int(C)  is  the  interior  of  convex  elastic  domain  and  we  see  from  (B.19)  that  in  the 

• • 

elastic  region,  X = 0 and  (f>int(X.)  = 0.  The  parameters  X cannot  evolve  unless  the  state 
of  generalized  forces  Y attains  a threshold  corresponding  to  the  boundary  of  the  convex 
domain , dlZ. 

The  above  normality  rule  states  that  the  flow  rule  may  be  written  in  the  form 


X = W(  Y) 


(B.22) 


for  a non-smooth  yield  surface,  or 


X = 7V0*(Y) 


(B.23) 
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if  it  is  smooth.  Here  7 is  a non-negative  scalar  which  is  required  to  satisfy  the  comple- 
mentary conditions 


7 > 0,  0*(Y)  < 0,  W(Y)  = 0 


(B.24) 


It  follows  that  the  set  of  generalized  stresses  satisfying  4>*(Y)  > 0 is  not  attainable. 
Now  if  4>*(Y)  < 0 then  the  behavior  is  elastic,  and  so  there  can  be  no  changes  in  the  inter- 
nal variables.  The  same  applies  if  Y lies  on  the  transformation  surface,  that  is  4>*{Y)  = 0, 
but  the  rate  of  change  of  Y is  such  that  the  internal  force  moves  into  the  elastic  region. 
This  is  referred  to  as  elastic  unloading.  Consider  also  the  situation  in  which  Y lies  on 
the  transformation  surface  and  remains  on  it;  this  is  referred  to  as  inelastic  loading,  and 
it  is  assumed  that  this  is  accompanied  by  a change  in  internal  variables.  In  summary,  we 
assume  that 


X = 0 if 


' df*{Y)  < 0 
< or 

k <f>*{ Y)  = 0 and 


30*(Y)  < 0; 


(B.25) 


X f 0 only  if  f*(Y)  = 0 and  df*( Y)  = 0.  (B.26) 

The  requirement  during  inelastic  loading  that 

dcj)*(  Y)  = 0 (B.27) 

is  known  as  the  consistency  condition. 

Dual  Properties. 

Mathematically  speaking,  the  dissipation  definition  (B.14)  is  equivalent  to  the  fol- 
lowing  expression  in  terms  of  X in  the  form 

P(X)  = sup  {Y**X}>0  (B.28) 

(Y*ec)  1 J 

From  the  use  of  the  definition  of  support  function22  (Figure  B.3)  in  which  Y is  the  point 
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Figure  B.3.  Support  function:  4>*{z)  = (x,  z ) is  the  support  function  of  the  convex  K at 

x. 


at  which  the  supremum  is  achieved,  the  term  in  the  bracket  is  the  rate  of  plastic  work 
or  dissipation.  We  may  notice  that  for  the  function  of  X,  such  that  V = f which  is 
differentiable  for  X f 0,  we  simply  obtain 


Y 


dV 

dk 


(or  similar  to  in  plasticity 


a cx 


df_ 

deP} 


(B.29) 


The  dissipation  potential,  T>(X),  if  associated  to  the  plastic  model,  is  a convex  function, 
positive  and  homogeneous  of  degree  one.23 

The  question  we  may  ask  now  is  the  following:  for  a given  Y,  try  to  find  X such 
that  the  Y might  be  associated  to  one  among  them  (Maugin  [1992,  p.53]).  Given  the 
properties  of  convexity  it  is  sufficient  to  introduce  the  notion  of  the  Legendre-Fenchel 
conjugate  V*  of  V by 


(B.30) 


22  inf  stands  for  infimum  or  minimum  and  sup  stands  for  supremum  or  maximum.  Simply,  ay  = inf  (.4) 
means  ay  is  a greatest  lower  bound  for  a set  A.  If  ai  = sup(T)  means  ai  is  a least  upper  bound  for  a set 
A. 

23  Then  <p  is  a positively  homogeneous  function  of  degree  one  of  the  flux  X,  that  is, 


</>(/fcX)  = kct>(k)  >0  if  k > 0 
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We  get  immediately 

9V*{Y)  = jx 


V*(Y*)  — V(Y)  > X*[Y*  - Y],  V Y” 


Consequently,  X is  a subgradient  of  V*  at  the  point  Y : 


X e dV*{Y)  = NC(Y) 


Therefore,  if  V*(Y)  is  continuously  differentiable  for  Y^O: 

dV* 


X = 


dY 


(B.31) 


(B.32) 


(B.33) 


The  properties  of  convexity  and  of  homogeneity  are  preserved  in  the  Legendre-Fenchel 
conjugate  V*(Y)  transformation  (B.30).  Eq.  (B.32)  is  the  dual  expression  to  (B.29) 
and  T>*{ Y)  may  be  said  to  be  the  dual  dissipation  pseudo-potential ,24  Eq.  (B.32)  or 
(B.33)  is  an  evolution  equation,  generally  nonlinear,  for  X.  More  detailed  discussions 
and  developments  concerning  Generalized  Standard  Materials  (GSM)  are  referred  to  as 
the  generalized  normality  hypothesis  and  to  Halphen  and  Nguyen  [1975]  and  Han  and 
Reddy  [1999], 


Definition  B.l.  Let  C C X,  a normed  vector  space.  The  interior  and  boundary  of  C 
are  denoted  respectively  by  int{C)  and  dC.  The  subset  C is  said  to  be  convex  if 

for  any  xx,x2  E C and  A G (0, 1)  , Xxx  + [1  - X\x2  G C (B.34) 

In  other  words,  the  subset  C is  convex  if  and  only  if  the  line  segment  joining  any  two 
points  of  C lies  entirely  in  C (Figure  B.5).  I 

Definition  B.2.  The  function  / is  said  to  be  convex  if 

f(Xxx  — [1  — A]x2)  < Xf(xx)  + [1  — A]/(x2)  Vx1(  x2  G X,  and  A £ (0, 1)  (B.35) 

The  function  / is  said  to  be  strictly  convex  if  the  strict  inequality  in  above  equation  holds 
whenever  xi  ^ x2  (Figure  B.6).  I 


24 


V*(Y)  is  none  other  than  the  yield  surface  fy  = 0 (Maugin  [1999,  p.  174]). 
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Figure  B.4.  Level  surface  0*(  er  , A ")  and  its  projection  onto  the  stress  space. 


Figure  B.5.  Convexity  of  a subset  C. 
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Figure  B. 6.  Convexity  of  a function  f(X). 
/(*) 


Figure  B.7.  Subgradient  of  a non-smooth,  convex  function. 


Definition  B.3.  Given  a convex  function  / on  a set  X,  for  any  x0  E X the  subdif- 
ferential df(x0)  of  / at  x0  is  subset  of  A'*  defined  by  (Reddy  and  Martin  [1994,  p.75]) 

df{x  0)  = jx*  6 X*  f(x)  - /(x0)  ><x*,x-x0>  y y E AT  j (B.36) 

A member  of  <9/(x0)  is  called  a subdifferent  of  / at  x0.  If  / is  differentiable  at  x0,  then 


df{x  0)  = {V/(x0)} 


(B.37) 


At  a corner  point  (x0,/(x0)),  the  subdifferential  <9/(x0)  is  the  set  of  the  slopes  of  all 
the  lines  lying  below  the  graph  of  / and  passing  through  the  point  (x0,  /(x0)).  This  is 
illustrated  in  Figure  B.7.  I 
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